Локальные и сплайновые аппроксимации в цифровой обработке геомагнитных наблюдений Текст научной статьи по специальности «Математика»

ЧЕБЫШЕВСКИЙ СБОРНИК

Том 19. Выпуск 4

УДК 550.380 DOI 10.22405/2226-8383-2018-19-4-26-42

Локальные и сплайновые аппроксимации в цифровой обработке

геомагнитных наблюдений 1

Гетманов Виктор Григорьевич — доктор технических наук, главный научный сотрудник, Геофизический центр FAH; главный научный сотрудник, Институт физики Земли им. О. Ю. Шмидта FAH. e-mail: v.getmanov@gcras.ru

Аннотация

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

Ключевые слова: локальные и сплайновые аппроксимации, модели, геомагнитные наблюдения, цифровая обработка.

Библиография: 22 названий. Для цитирования:

В. Г. Гетманов. Локальные и сплайновые аппроксимации в цифровой обработке геомагнитных наблюдений // Чебышевский сборник, 2018, т. 19, вып. 4, с. 26-42.

1Статья написана в рамках Государственного задания ГЦ РАН на 2018 г.

CHEBYSHEVSKII SBORNIK Vol. 19. No. 4

UDC 550.380 DOI 10.22405/2226-8383-2018-19-4-26-42

Local and spline approximations in digital processing of geomagnet-ic observations

Getmanov Viktor Grigorievich — D.Sc., Principal research scientist, Geophysical Center RAS; Principal research scientist, Schmidt Institute of Physics of the Earth RAS. e-mail: v.getmanov@gcras.ru

Abstract

Methods of local and spline approximations for digital processing of geomagnetic observations are proposed for consideration. Algorithms for calculating piecewise-linear, sinusoidal and pol-ynomial local approximation models have been developed. An algorithm for calculating the spline approximation model is developed. The generated mathematical apparatus is focused on solving problems of parameter estimation, filtering and spectral analysis for geomagnetic observations.

Keywords: local and spline approximations, models, geomagnetic observations, digital processing.

Bibliography: 22 titles. For citation:

V. G. Getmanov, 2018, "Local and spline approximations in digital processing of geomagnet-ic observations" , Chebyshevskii sbornik, vol. 19, no. 4, pp. 26-42.

1. Введение

В настоящее время обсерваторские геомагнитные наблюдения (ГМН) подвергаются, в основном, первичной обработке для решений задач рутинных статистических вычислений или редактирования данных; одновременно, действующая научно- техническая практика выдвигает вполне очевидные требования к повышению качества ГМН, расширению математического арсенала средств вторичной цифровой обработки и построению глобальных геомагнитных моделей. Создание новых методов цифровой обработки является актуальным направлением геофизики и геоинформатики [1-4]. Существенный вклад в развитие указанного направления внёс академик РАН Гвишиани А.Д. [5-8].

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

Предлагаемые в статье методы и алгоритмы локальных и сплайновых аппроксимаций для цифровой обработки в предметной области ГМН, предназначаются для задач оценивания многопараметрических моделей, фильтрации и сглаживания сложных зависимостей и реализации обобщённого спектрального анализа. Математический аппарат, рассмотренный здесь в обзорном плане, базируется на результатах работы [9].

Представленные здесь методы и алгоритмы нацелены на задачи цифровой обработки ГМН, которые традиционными подходами решаются не в полной мере эффективно. Так, содержащиеся в распространённых программных пакетах стандартные цифровые фильтры плохо работают на ограниченных временных интервалах; спектральный анализ на основе дискретного преобразования Фурье реализуется для нестационарных ГМН с большими погрешностями и

т д

2. Модельные представления для геомагнитных наблюдений

2.1. Обычно, компоненты вектора напряжённости геомагнитного поля измеряются в земной системе координат, для которой ось координат ОХ направлена на северный полюс, ось ОУ- по параллели, ось 02- по радиусу к центру Земли. В практике анализа характеристик геомагнитного поля достаточно часто координаты и модуль вектора напряжённости НХ(Ъ) = Н\(¿), Ну(¿) = Н2(Ь), Нх(¿) = Н3(Ъ) и Н(Ь) = ЯоСО представляются исходными функциями времени заданного вида.

Как правило, характеристики геомагнитного поля определяются на заданном временном интервале (tо,tf). В ряде случаев, зависимости от времени для координат и модуля вектора напряжённости геомагнитного поля представляются низкочастотными функциями времени. Для некоторых задач цифровой обработки определяется множество- класс функций, к которым должны принадлежать исходные функции времени. Наиболее подходящим, для этой цели является множество С(tо,tj); функциональные характеристики геомагнитного поля, которые принадлежат к С(tо,tj), определены на (tо,tf), вместе с производными вплоть до Ьо-ото порядка. Возможны наложения ограничений на производные для исходных функций времени.

2.2. ГМН координат Уп(Тг),п = 1, 2, 3 производятся векторными магнитометрами, наблюдения модуля Уо(Тг) производятся скалярными магнитометрами, Т - интервал дискретизации. Большая часть обсерваторий системы 1.\ТК1>.\1.\(!.\Т.Т. производящих ГМН с минутной дискретизацией- (Т = 60А), оснащены векторными индукционными магнитометрами УМ300 и скалярными магнитометрами БМЭООЫ фирмы (¡КОМЛС (Франция). Обсерватории, производящие ГМН с секундной дискретизацией (Т = 1А), большей частью оснащаются векторными

магнитометрами VM391 и скалярными магнитометрами фирмы GEOMETRICS G856. ГМН с минутной дискретизацией в стандарте INTERMAGNET доступны на сайте [10]; наблюдения с секундной дискретизацией - на сайте [11]. Бывают ГМН с неравномерной дискретизацией, дискретизацией по часам, суткам и т.д.

2.3. Пусть ГМН Yn(ti) в общем случае для неравномерной дискретизации to,t^t2,... описываются следующими модельными соотношениями

Yn (ti) = Hn(ti) + Wns(ti) + Wnr (ti),

где wns(ti) -систематические низкочастотные приборные погрешности, wnr(ti) -случайные погрешности; имеет место соотношение wnr(U)= wnr1(ti)+wnr2(ti) , где wnr1(ti) -случайные приборные погрешности, wnr2(U) -случайные погрешности, имеющие техногенное или природное происхождение. Модель наблюдений для функции модуля имеет вид

Yo(ti) = Ho(ti)+ wo(ti),Ho(ti) = ^ H2(ti) + H22(ti)+ Hl(ti).

Погрешности наблюдений wnri (ti),Wnr2(ti),wo(ti), обычно, моделируются в виде некоррелированных случайных чисел с нулевыми математическими ожиданиями и заданными дисперсиями.

Приведём величины погрешностей векторных и скалярных магнитометров. Для минутных наблюдений предельная погрешность магнетометра VM300 составляет % 5пТ, погрешность SM900R- ~ 1пТ. Для секундных наблюдений для векторного магнитометра VM391 предельная погрешность составляете (1 ^ 2)пТ, для магнитометра GEOMETRICS G856 - составляет величину ~ (0.5 ^ 0.7)пТ.

Отметим, что обсерваторские ГМН бывают двух видов: 1. первичные preliminary data, которые не были подвергнуты предварительной цифровой обработке; 2. definitive data - данные, подвергнутые цифровой обработке, с устранёнными систематическими погрешностями и уменьшенными случайными погрешностямии.

3. Методы локальных и сплайновых аппроксимаций

3.1. Рассмотрим ij - неравномерные моменты дискретизации, г = 0,1,...,N — 1 на локальном интервале (to,tN-1) и Y(ti) - произведённые ГМН. Поставим этим наблюдениям в соответствие локальную модель Нм(c,ti) в виде известной нелинейной функции, зависящей от вектора параметров ст = (с1,..., ст), в общем случае А eCi . Пусть-ш(^) - погрешности наблюдений, являются случайными нормальные числа с нулевым математическим ожиданием и дисперсией а2. Связь наблюдений, модели и погрешностей представим соотношением

Y (ti ) = Нм (c,U)+ w(ti).

Сформируем локальный функционал Si(Y,c), реализуем его условную минимизацию, найдём с°

N-1

Si(Y,c)= V(Y(ti) — Нм(c,ti))2,A° = arg{min Si(Y,c)}. (1)

Последовательность H^L(ti) = Нм(c°,U), i = 0,1,...,N — 1 для c° из примем в качестве результата реализации метода локальной аппроксимации.

3.2. Рассмотрим исходный временной интервал (to,tNf- моменты времени наблюдений, Y(ti),i = 0,1,...,Nf — 1 - ГМН на данном интервале. На исходном интервале зададим сплайновые узлы Ti, Т2,..., тп-1 и сформир уем п сплайновых интервалов -

(t0 < t ^ т1), (т1 < t ^ т2),..., (тп-1 < t < Nf — 1). Введём интервальные дифференцируемые функции f (а к, t), которые равны нулю вне сплайновых интервалов, ак, к = 1,...,п -векторы параметров, ат = (af, ...,&П)• Зададим ограничивающее множество

Ä2 = {« : (f(s)(ак,тк) — f {s)(ak+i, Тк) = 0,к = 1, ...,п — 1, s = 0,1)}. Определим сумму интервальных функций

п

f (a,t)= ^ f (ак,t). к=1

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

Сформируем сплайновый функционал S2 (Y, а), решим задачу его условной минимизации для отыскания вектора оптимальных параметров а°, обеспечим нахождение оптимальной сплайновой аппроксимационной модели

N-1

S2(Y,a) = V (Y(ti) — f (a,ti))2,a° = arg{ min S2(Y,a)}. (2)

7=o aeA

Последовательность H^s(ti) = f (a°, ti), i = 0,1,..., Nf — 1 для a° из (2) примем в качестве результата реализации метода сплайновой аппроксимации.

В [12] помещены более подробные материалы по постановкам сформулированных здесь локальным и сплайновым аппроксимациям; работа [13] содержат обобщения, касающиеся технологии локальной аппроксимации. Методы п.п. 2.1, 2.2 базируются, в общем случае, на применении нелинейного программирования и представляют собой обобщения подходов традиционных регрессионных моделей. Их особенности состоят в построении аппроксимационных моделей на малых локальных интервалах, учёте связей моделей на соседних локальных интервалах, использовании локальных моделей упрощенного вида с малым числом параметров и возможности аппроксимации структурно сложных и существенно нестационарных ГМН.

4. Построения локальных аппроксимационных моделей

Здесь представим варианты построений локальных аппроксимационных моделей, используемых при реализации цифровой обработки ГМН.

4.1. Рассмотрим построение локальных аппроксимационных моделей, линейных по части параметров. Положим, что задан локальный интервал с точками г = 0,1, ...,М — 1 и К(Тг) -ГМН. Для указанных локальных моделей, вектор параметров ст = ([5Т,шт) является блочным и состоит из векторов @ и ш размерности тю, ш2о, которые в модель входят линейно и нелинейно. Запишем соотношения для модели

тю

Нм(с, Тг) = ^ РГ1 <рГ1 (ш,П)= ¡3Тч>(ш,Тг),

Г\ = 1

где ^(ш, Тг) - векторная базисная функция, зависящая от вектора параметров ш

ф, Тг)т = (<р1(ш, Тг), ч>2(ш, Тг),..., ^(и, Тг)).

Как правило, на вектор @ не накладываются ограничения; координаты вектора ш являются ограниченными - ш € О,о.

Введём функционал 631 (У, Р, ш) и запишем основную задачу минимизации по ш, которую будем решать в два этапа

N-1

N-1

S3i(Y,ß,u)=J2(Y(Тг) - Нм(ß,u,Ti))2 = (тi) - ß<P(u,Ti))

i=0

i=0

(в°,u0) = arg{ min S31(Y,f3,u)}.

На первом этапе фиксируем допустимое ш = const и будем находить частично оптимальные линейные параметры, зависящие от ш

¡3°(u)=arg{ mm S3i(Y,p,u)}.

Р,ш=с onst

Сформируем вектор наблюдений У, матрицу плана X(ш) и вектор параметров ft ,Х (и) =

Y =

Y(Т • 0) Y(Т • 1)

Y(Т(N - 1))

<Pi(w,T • 0) • 0)

Vi(u,T • 1) V2(u,T • 1)

(fimw(ш,Т • 0) (ш,Т • 1)

Vi(u,T(N - 1)) v2(u,T(N - 1)) . (u,T(N - 1))

ß(u) =

ßi(u) ß*(u)

ßmw

Воспользуемся [14], запишем систему линейных уравнений с использованием введённых обозначений А(ш)Р°(ш) = Ь(ш), где Ъ(ш) =ХТА(ш) =ХТ(ш)Х(ш). Решение данной системы позволит найти частично оптимальный вектор @°(ш).

На втором этапе подставим @°(ш) в выражение для функционала , получим выражение для частичной остаточной суммы, зависящей от ш

Sa,3i(Y,ß°(ш),ш) = YTУ - b1 '(u)ß°(ш).

иТ/

(3)

Найдём оптимальное значение вектора ш° па основе минимизации частичной остаточной суммы Бо,з1(У, Р°(ш),ш) (3), численное значение которой можно вычислить для любого допустимого ш

2

= arg{ min S0,3i(Y, ß°(w), w)}.

шбПо

Минимизацию реализуем на основе поисковой процедуры для ш £ Q,q. Оптимальный вектор

т т т

ß° найдём из решения системы A(u°)ß° = Ь(ш°), получим с° = (ß° ,ш° )и выражение для оптимальной модели Нм(с°,Тг)

4.2. Рассмотрим построение аппроксимационных моделей на системе локальных интервалов. Произведённые ГМН представим в виде У (Тг), г = 0,1,...,Nf — 1. Исходный интервал разобьём на m локальных интервалов по Жточек, N(j — 1) + 1 ^ i ^ Nj, mN = Nf. Введём локальные модели Нм (cj, Тг), зависящие от параметров Cj, j = 1,...,m , Нм (cj ,Ti) = 0 для г < N(j — 1), г > Nj . Зададим локальные функционалы и основной функционал

Nj-l m Nj

S32 (Cj, н )= £ (Hj (Ti) — Нм (Cj .Ti))2, S32 (с, н ) = £ £ (H (Ti) — Нм (cj, Ti))2

i=N(j-l) j=l i=N(j-l)+l

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

с° = arg{ min S2(Cj, H)}, j = 1,..., т.

Cj € Co

Оптимальная аппроксимационная модель на системе локальных интервалов представляется последовательностью кусочно- непрерывных функций и имеет вид

т

Н°м(Ti) = ^ Нм(cj,Ti), i = 0,1,..., Nf — 1. j=i

4.3. Рассмотрим вычисление скользящих локальных моделей и процедуру их взвешенного усреднения, которые реализуем в два этапа [15].

На первом этапе для исходного интервала с точками г = 0,1,..., Nf — 1 введём N -точечные скользящие локальные интервалы с граничными точками Nij,N2j, Л^-шаг скольжения, то-число скользящих интервалов

Nij = Nd(j — 1),N2j = Nij +N — 1,j = 1,..., то.

Для упрощения потребуем выполнения соотношений кратности mN = Nf и Ndmd = N, т, md-целые числа, тогда то = md(m — 1) + 1. Для скользящего локадьного интервала с номером j найдём модели Н( cj ,T г), j = 1,..., то, где cj - оптимальные параметры скользящих локальных моделей.

На втором этапе введём единичные функции Е-(Tг), равные нулю вне локальных скользящих интервалов, вычислим их сумму Е(Tг) и функцию последовательности весовых коэффициентов R(i)

Ej(Ti) = 1, Nij < i < N2j,Ej(Ti) = 0, 0 < i < Nij, N2j <i < Nf — 1,

то

E(Ti) = ^ Ej(Ti),R(Ti) = 1/E(Ti),i = 0,1,..., Nf — 1. (4)

=i

Произведём взвешенное усреднение с использованием (4) для суммы последовательности скользящих локальных моделей.

то

Н°м (T i) = R(i) £Нм (С j, T i). =i

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

Локальные кусочно- линейные аппроксимационные модели будем считать базисными для задач цифровой обработки ГМН. Рассмотрим алгоритм построения указанных моделей.

По- прежнему, У(Тг) - произведённые ГМН на локальном интервале, г = 0,1,..., N — 1. Локальную модель запишем следующим образом: Нм(с, Тг)= с1 + с2Т1, где ст = (с1, с2) - вектор параметров. Базисные функции представим векторами размерности ( N, 1) <р1(Тг)Т = (1,..., 1), !^2(Тг)т = (Т■ 0,Т■ 1, ...,Т(М — 1)); модель примет вид Нм(с, Тг)= сТ^(Т1). Погрешности

примем случайными нормальными числами с нулевым математическим ожиданием и дисперсией а2. ГМН, модель и погрешности связаны соотношением

У (Т г) = стр(П)+-ш(Т г).

Ставится задача по наблюдениям У ("г), г = 0,1,...,И — 1 найти оптимальные оценки параметров с1, с2 и оптимальную локальную кусочно- линейную модель. Сформируем локальный функционал и выражение для задачи минимизации

N-1

Я*!(У, с) = V (У(Тг) — ст^(Тг))2,с° = с)}.

Введём матрично- векторные обозначения для записи необходимых выкладок для задачи минимизации Ут = (У(Т ■ 0),У(Т ■ 1),...,У(Т(И — 1))) X = (^(Тг), <р2(Т1)), А = ХтХ,

Ь = X тн.

А А =

А,

N-1 N-1

агз = ^ (Тг)(р8(Тг)г, в = 1,2, Ъг = ^ (Тг)Н(Тг),

г=0 г=0

N-1 N-1 N-1 N-1

ап = Иа12 = Т ^ га21 = а^2 = Т2 ^ г%1 = ^ Н(Тг), Ь2 = ^ Н(Тг)Тг.

г=0 г=0 г=0 г=0

Вычислим решение линейной системы (с1, с2)~, запишем выражение для оптимальной локальной аппроксимационных кусочно- линейной модели

апс°1 +аис2 = Ь1,а21С°1 + а22с2 = Ь2 ,Н°М (с° ,Тг) = с\ + с%П.

6. Алгоритмы вычисления локальных синусоидальных и полигармонических аппроксимационных моделей

6.1. Рассмотрим алгоритм вычисления локальных кусочно синусоидальных аппроксимационных моделей [16]. Пусть У(Тг) г = 0,1,..., N — 1 -ГМН. Примем в качестве локальной модели функцию Нм(с, "г) = асоншТг + Ь8\ишТг в предположении, что функция геомагнитного поля на локальном интервале является одночастотной. Амплитуды а, Ь входят в модель линейно, частота ш - нелинейно; положим, что амплитуды не подвержены ограничениям, частота ш ограничена диапазоном О0 = (ш : штщ ^ ш ^ штах), со значениями шт\пж штах. Допустим, что

погрешности наблюдений w(Тг)-случайные некоррелированные нормальные числа с нулевым

а2

У (Тг ) = Нм (с, "г)+w(Тг).

а ,

частоты ш° для локальной кусочно- синусоидальной аппроксимационой модели. Функционал 351(У,а,Ъ,ш) сформируем следующим образом

N-1

8ы(У,а,Ь,ш) = ^ (У (Т г) — а еовшТ1 — ЬвтшТ I)2. =о

Оценки параметров модели найдём путём минимизации по а,Ь, ш € Оо

(а°,Ь°,ш°) = arg{ min S(Y,a,b,w)}.

a,b,oj£Qo

Запишем модель в виде Нм(с,"i) = ßT<ß(w,Ti), ßT = (а, b). Базисные функции ^i(Ti), (р2(Гх)представим векторами размерности (N, 1): ipi(Ti)T = (cos(uT • 0),cos(uT • 1),... ..., cos(uT (N-1)), ^2(Ti)T = (sin(uT •O), sin(uT •l),..., sin(uT(N—1)). Воспользуемся векторпо-матричными обозначениями YT = (Y (T • 0),Y (T • 1), ...,Y (T (N — 1))) X = (p1(Ti), <p2(Ti)), A = XTX, b = XTH. Решение задачи минимизации осуществим в два этапа.

На первом этапе для фиксированной частоте ш £ Qo найдём решение линейной системы А(ш)А°(ш) = Ь(ш) . Выпишем коэффициенты для А(ш),Ь(ш)

N-1 N-1

агs(ш) = ^ (и, Ti)<p8(u, Ti)г, s = 1,2, br(ш) = ^ (и, Тг)Н(Ti),

г=0 г=0

N-1 N-1 N-1

а11(ш) = ^ cos2 wTi, 012(w) = ^ cos wTi sin wTi,a21(w) = а12(ш), а22(ш) = ^ sin2 wTi, =o =o =o

N-1 N-1

h(u) = ^ Y (Ti) cos uTi, Ъ2(ш) = ^ Y (Ti) sin wTi.

i=o i=o

Оптимальные амплитуды найдём из решения линейной системы уравнений

а°(ш)ап(ш) + Ь°(ш)а,12(ш) = h(w), а°(ш)а,21(ш) + Ь°(ш)а,22(ш) = Ь2(ш),

а°(ш) = (Ъ\(ш)а22(ш) — Ь2(ш)а\2(ш))/А(ш),Ь°(ш) = (Ь2(ш)ац(ш) — Ь\(ш)а12(ш))/А(ш)

где А(ш) = а11(ш)а22(ш) — а'22(ш) - определитель.

На втором этапе сформируем частичную остаточную сумму

N-1

SoMY, ш)= £ Y2(Ti) — а°(ш)Ьг(ш) — Ь°(ш)Ь(ш).

г=0

Нахождение оптимальной оценки частоты ш° реализуем путём минимизации So;5i(Y, ш) на основе подпоиска по частоте ш в диапазоне Qo

ш° = arg{ min So;5i(Y,w)}.

шбПо

Подставим в формулы а°(ш), Ь°(ш) значение ш = ш°, найдём искомые оптимальные оценки амплитуд а° = а°(ш°), Ь° = Ь°(ш°) и запишем выражение для оптимальной модели Нм(с°,Тг) = а° cos w°Ti + b° sinu°Ti.

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

6.2. Рассмотрим алгоритм вычисления локальных полигармонических аппроксимационных моделей для локального интервала с точками i = 0,1,..., N — 1 с ГМН Y(Ti) для исходной многочастотной функции геомагнитного поля. Аппроксимационную модель в предположении, что исходная функция геомагнитного поля на локальном интервале является многочастотной, представим в виде соотношения

L

Нм(с, Ti) = cos wiTi + bi sin wiTi),

=i

где вектор ст = (ат, Ьт, шт) = (а1,а2,..., аь, Ь1, Ь2,..., Ьь, (1,(2, ...,шь) имеет размерность (3Ь, 1) и состоит из векторов амплитудных и частотных параметров. На линейные амплитудные параметры ограничения не накладываются: —те ^ а ^ те, —те ^ Ь1 ^ те, 1 = 1, 2,..., Ь. Примем, что исходная многочастотная функция геомагнитного поля реализована в ограниченной полосе частот с границами штт,штэ,х~, для векторов нелинейных частотных параметров с учётом их упорядочивания, допустимое множество представим следующим образом

= {(^1, ...,Шь) : ¿т1п ^ ¿1 <¿2 < ... <Шь ^ ¿тах},( € ^ Запишем функционал

N-1 N-1

\2

S52(Y,a,b,w) = (Ti) -Нм (c,Ti))2 (Ti) (ai COSWl Ti + bl s[nWl Ti))2-

i=0 i=0 1=1

Вычислим оптимальные параметры аппроксимационной модели с помощью решения задачи минимизации введённого функционала

(a°,b°,w°) = arg{ min S52(Y,a,b,w)}.

а,Ь,ш£.£1 q

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

L

л; jcos =1

ai coswi Ti + bl sinwiTi) = ßTp(w,Ti),i = 0,1,...,N - 1,

где вектор линейных параметров ßT = (ai,..., aL, bi,..., bL) имеет размерность (2L, 1) ^(w, Ti) - векторная базисная полигармоническая функция

<f(w, Ti)T = (cos(w1Ti),..., cos(wLTi), sin(w1Ti),...., sin(wLTi)).

Пусть вектор ГМН Y имеет размерность (N, 1). С использованием <ß(w,Tг) сформируем матрицу плана %(w) размерности (N, 2L)

YT = (Y(T ■ 0),Y(T ■ 1),..., Y(T(N - 1)), %%(w) = \<p(w, T ■ 0),..., p(w, T(N - 1))|T .

Тогда функционал запишем в векторно - матричном виде:

S52(Y, ß, w) = (Y - X(w)ß)T(Y - X(w)ß),

На первом этапе минимизации для S(Y, ß, w) фиксируем нелинейные параметры w = const, w E &L, находим частично оптимальные линейные параметры ß°(w)

ß°(w)=arg{ min S52(Y,ß,w)},

ß,Lü=c an st

которые вычислим на основе решения системы линейных уравнений

A(w) = XT(w)X(w), b(w) = XT(w)Y, ß°(w) = A(w)-1 b(w), (5)

На основе ß°(w) сформируем функционал частичной остаточной суммы, зависящий только w

SoMY,u) = S(Y,ß°(w),w) = YTY -ß°T(w)b(w),

На втором этапе аппроксимации для So (Н, ш) осуществим минимизацию по w £ Qописание поискового алгоритма содержится в [18] Находим оценки частотных параметров ш° и, после подстановки ш° в (5), вычислим оценки амплитудных параметров ß° для локального интервала с точками г = 0,1, ...,N — 1 и запишем выражение для оптимальной модели

L

ш° = arg{ min So (Y,w)},ß° = ß° (w°) = A(w°)-1 b(w°), Н(с°,Тг) = J^(af cosw^Ti + bf sinw^Ti).

Публикация [19] содержит материалы, касающиеся приложений вычислений локальных апппроксимационных полигармоничексих моделей для реализации обобщённого спектрального анализа ГМН.

7. Алгоритм вычисления локальных аппроксимационных моделей для цифровой обработки ГМН от векторных и скалярных магнитометров

7.1. Положим, что Y\(Ti), Y2(Ti),Ys(Ti)rn Yq(Ti) - ГМН координат и модуля вектора напряжённости геомагнитного поля, произведённые векторными и скалярными магнитометрами на локальном интервале, г = 0,1, ...,N — 1. Определим для координат вектора напряженности локальные модели на основе нелинейных известного вида функций Hi(ci,Tг), Н2(c2,Tг), Нз(C3,Ti), где Ап -векторы параметров моделей, п = 1, 2, 3. Эти модели могут быть, в частности, кусочно -постоянными, кусочно- линейными и т.д. Будем считать, что погрешности наблюдений относятся к definitive data -в них отсутствуют систематические погрешности (уменьшены в максимально возможной степени) Представим wn(Ti),n = 0,..., 3 в виде случайных некоррелированных нормальных чисел с нулевыми математическими ожиданиями и дисперсиями а2. Свяжем наблюдения, модели и погрешности линейными зависимостями

Yn(Ti) = HMn( Cn,T i) + Wn(Ti), п = 0,..., 3.

Ставится задача по definitive data Yn(Ti),n = 0,1, 2, 3 нахождения оценок локальных моделей H°Mn(Ti) = Нмп(cl,Ti).

Основываясь на определенных выше наблюдениях и моделях, локальный функционал ^6i(Y, с), определяющий меру близости локальных наблюдений и моделей, запишем в виде суммы четырёх функционалов

N-1

So,ei(Yo, с) (Ti)

i=0 N-1

£нМп(с n ,T i))2

n=1

Sn,6i(Yn, c) =Y,(Yn(Ti) -YMn(Cn,Ti))2,n = 1, 2, 3

i=0

в61 (У, с) = БомУо, С1, с2, сз) + ,61 (Yп, Сп),

п=1

где ст = (с^, с2 , ), Yт = (ут , Y2[, Y3f). Нахождение оптимальных оценок параметров моде-

п

задачи оптимизации локального функционала

с° = &щ({тт361(У, с)},с°Т = (, <£ , с0^).

Очевидно, что на локальном интервале исходные функции координат должны быть почти постоянными. Определим, в этом случае, локальные модели в виде кусочно- постоянных функций Нмп(сп, Тг) = Ап, п = 1, 2, 3 и выражение для локальных функционалов

0 N-1 N-1

Бегр, с) = £(£ (¥п(Тг) - <п)2 + £ (Уо(Тг) с2 + с2 + с2)2. (6)

В=1 г=0 г=0

Ап

экстремума в виде системы из трёх нелинейных алгебраических уравнений

= Т(УпТ) - 2»с + АпШ+?) = = 1■ 2• 3 <7>

и Ьп г=0 УС1 + С2 + С3

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

л •¡г^N-1 у (Т ) N-1

< = 7^(1 + , Ьг=оУо(Т V == у Уп(П),п = 1, 2, 3.

п 2» ..—.М 1 ______________„ ЛГ !__._..„ „ ЛГ 1 __._..„ ' у '

2N" ' \/(Zr=-0 Y1(Tr))2 + (Zr=-0 У2(Тг))2 + (Zr=-0 Ys(Tr))2 г=0

Оценки локальных аппроксимационных моделей представим в виде H'^Ln(Ti) = с п = 1, 2, 3.

Описание результатов и приложений п.п. 6.1. помещено в публикации [20].

7.2. Положим, как и в п.п. 6.1, что Yn(Ti), п = 0,..., 3 - ГМН, произведённые векторными и скалярными магнитометрами на локальном интервале, i = 0,1,..., N — 1. В данном случае погрешности наблюдений wn(Ti),n = 0,..., 3 отнесём к preliminary data -первичным данным, которые не были подвергнуты предварительной компьютерной обработке и в которых присутствуют, как систематические, так и случайные погрешности. Обычно, для векторных магнитометров на локальных интервалах медленно меняющимся систематическим погрешностям поставим в соответствие модельные кусочно- постоянные функции - Сп, п = 1, 2, 3; в общем случае, считаем, что с £ С. Для скалярных магнитометров доиустим, что Со = 0. Погрешности wn(Ti) представим в виде некоррелированных нормально распределённых случайных чисел с нулевыми математическими ожиданиями и дисперсиями аП, п = 0,1, 2, 3, причёма2 <<(тП

п = 1, 2, 3 мостями

Yn(Ti) = Hn(Ti) + Cn + Wn(Ti), n = 0,..., 3.

Ставится задача no preliminary data Yn(Ti),n = 0,1, 2, 3 нахождения оценок модельных систематических погрешностей с^, с2,с% .

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

N-I

S62(Y, с) = £ (Y2(Ti) — (Yi(Ti) — ci)2 — (Y2(Ti) — C2)2 — (Y(Ti) — c3)2)2.

i=0

Оценивание систематических погрешностей на локальном интервале реализуем на основе следующей задачи условной минимизации

т

с° = а^{шт вб2(У, с)}, с° = (с°1, с°2, с°). се с

8. Алгоритм вычисления сплайновых аппроксимационных моделей для ГМН

Здесь ставится задача по произведённым ГМН У(Тг) вычислить, в соответствие с 1.2., оптимальную сплайновую аппроксимационную модель (Тг) , г = 0,1,...,Nf — 1. Для

У( Т )

локальных моделей из п.п.3.2. Решение осуществим на основе задачи минимизации квадратичного функционала с линейными равенствами- ограничениями. Подходы к её решению в деталях описаны в [21, 22].

Введём на исходном интервале п сплайновых интервалов по N точек, N4 = Nf. г-ые точки к-го сплайнового интервала подчиняются неравенствам N (к — 1) ^ г ^ (^к — 1), к = 1, ...,п. Точки Nк , к = 1, ...,п — 1 примем в качестве стыковочных.

ций Д (1) =( /о,к (^, /1,к &),..., /ь,к &))Т размерности Ь, которые могут быть, в частном случае,

к

векторы весовых коэффициентов а^ = (ао,к,а1,к,...,аь,к), интервальную функцию для к-го интервала представим в виде /Т(1 )ак• Дискретизируем интервальные функции, просуммируем их и обозначим сумму как /Зр(а, Тг) , аТ = (аТ, аТ,..., аТ)

и(а, Тг) = ^ /к (Тг )ак, 0 < i < Nf — 1. к=1

Введём функционал Бт (У, а) , который является квадратичным

п —к-1

Бт(У,а) = £ £ (У (Тг) —/Т (Тг )ак )2,

( 1 ( Т <■) — ^ ( Т <■ )а^ 2

к=1 г=— (к-1)

Бт(У, а) =Т,П=1[Т,2кщк-1) (У 2Т) — 2 Е=-А-1) У (Тг (Тг )ак + Е—=к-к-1) (Тг (Тг)

Введём матрицы Qк, к = 1,..., п, которые, исходя го , имеют размерность ((Ь + 1), (Ь + 1))

- к-1

Qк = Е )]т (Т1).

г=—(к-1)

Для вектор-столбцов Рк, к = 1,..., п размерноети (Ь + 1) вычислим элементы Рк,1

-к-1 -к-1 Рк = Е У(Т1 )Г(Тг),Рк,1 = Е У (Т1) МТг),1 = 1,...,Ь,Рк = (Рк„Рк,1 ,...,Рк,ь)

г=—( к-1) г=—(к-1)

Рк

п —к-1 п п

Бт(У,а) = Е Е У2(Т{) + ^®ТQкак — 2 ^Р^к. к=1 —( к-1) к=1 к=1

Т

Образуем блочно-диагональную матрицу ф размерности ((Ь + 1)п, (Ь + 1)п) и блочный вектор РТ = (РТ, Р2^,..., Рп ) размерности п(Ь + 1), представим Бч(У, а) в векторно-матричной форме

п N^1

Бт(У, а) = ^ ^ У2(Тг) + аТфа - 2Рта. (8)

к=1 г=N (к-1)

Гладкость для , по нулевой и первой производным в стыковочных точках обеспечим с помощью линейных по а равенств кок(а)= 0 кк(а)= 0, к = 1, ...,п - 1

Ьок(а) = ¡Т(Т»к)ак - /Т+1(Т»к)ак+1 = 0,1цк(а) = ¡к1)Т (Т»к)ак - ¡к+1(Т»к)ак+1 = 0. (9)

Регулирование нулевыми и первыми производными (положением и наклоном) на концах интервалов наблюдения реализуем с помощью задания значений нулевых и первых производных на концах в точках Т ■ 0 и Т»п осуществим с помощью вектора д = (д 1, д2,до, д4)Т и линейных по а равенств к2з(а) - д = 0, з = 1,..., 4

Ыа) = аТ¡1 (Т ■ 0) - дх = 0, к22(а) = аТ/(1)(Т ■ 0) - 92 = 0,

к2о(а) = аТ1п(ТМп) - до = 0, Ъы(а) = ^¡^(ТХп) - д4 = 0.

(10)

Для (9) введём матричные условия-равенства. Блочные матрицы Оо,Б1 с размерностями ((п - 1), (Ь + 1)п) сформируем из векторов ¿оо = ¡Т(ТМк), (о1 = -¡Т+1(Т»к),

(1)т (1)т Л1о = ¡к ) (ТN к) , (11 = (ТN к). Для сформируем матрицу Б2 размерности (4, (Ь + 1)п),

образуем векторы (21 = ¡Т(Т ■ 0^ ((22 = Ц4' (Т ■ 0^ (20 = ¡Т'(Т»п), (24 = Ц4' (Т»п)

(1)

Т!

Л1)1

Бо

оо (о1 0 0 1о (11 0 0

0 ¿оо ¿о1 0 0 ,О1 = 0 ¿1о (11 0 0

0 оо (1о1 0 0 1о (11 0

0 0 (оо (о1 0 0 (1о ¿11

П2 =

(21 0 (22 0

0 ■■■

0 ...

. . . 0

. . . 0

0 (20

0 ¿24

(11)

Линейные равенства (9), (10) с использованием (11) запишем в матричном виде

Боа = 0, О1а = 0, Б2а - д = 0.

(12)

Решение задачи оптимизации Бт(У, а) (8) с учётом условий (12) для вектора а осуществим с помощью Хо, Х1, Х2 и Ь(У, а, X)- множителей и функции Лагранжа

п N1^-1

Ь(У,а,Х) = ^ ^ У2(Т1)+аТфа - 2РТа + ХТОоа + ХТО1а + ХТОа - д). (13)

к=1 г=N (к-1)

Ь( У, а, Х)

а Хо Х1 Х2

2фа - 2Р + БТХо + оТХ1 + БТХ2 = 0, Боа = 0, Бха = 0, Б2а -д = 0. (14)

Найдём решение линейной системы (14) а° и с применением вычислим оптимальную сплай-новую аппроксимационную модель H^s(Ti) =f°p(a°,Ti), i = 0,1, ...,Nf — 1.

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

9. Заключение

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

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

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

СПИСОК ЦИТИРОВАННОЙ ЛИТЕРАТУРЫ

1. Гвишиани А.Д., Лукьянова Р.Ю. Геоинформатика и наблюдения магнитного поля Земли: российский сегмент // Физика Земли. 2015.№2. С.3-20. doi:10.7868/S0002333715020040.

2. Гетманов В.Г. Цифровая обработка нестационарных колебательных сигналов на основе локальных и сплайновых моделей. М.: Изд-во НИЯУ МИФИ. 2011.298с.

3. Гетманов В.Г., Сидоров Р.В., Дабагян P.A. Метод фильтрации сигналов с использованием локальных моделей и функций взвешенного усреднения.//Измерительная техника. 2015. Ш. С.52-57.

4. Гетманов В.Г. О частотном подпоиске в задаче оценивания параметров кусочно- синусоидальных функций // Автометрия. 1992. № 2. С. 93-98.

5. Гетманов В.Г., Дабагян P.A., Сидоров Р.В. Исследование характеристик геомагнитных пульсаций методом локальных аппроксимаций // Геомагнетизм и аэрономия. 2016. № 2. Т.56. С.209-216.

6. Гетманов В.Г. Об алгоритме поиска по частоте в задаче оценивания параметров моделей полигармонических сигналов //Автометрия.-2009.№3.С.83-89.

7. Гетманов В.Г., Борзунов Г.И. Алгоритм параллельных вычислений для задачи спектрально- временного анализа на базисных полигармонических функциях // Информационные технологии. 2015. Т.21. №6. С.456-463.

8. Гетманов В.Г. Нелинейная фильтрация наблюдений системы векторного и скалярного магнитометров // Измерительная техника. 2013. №6. С.51-55.

9. Гетманов В.Г. Методы вычисления аппроксимационных сплайновых функций для задач цифровой фильтрации // Теория и системы управления. 2016. №5. С.45-53.

10. Гетманов В.Г. Алгоритмы вычисления аппроксимационных сплайновых функций с учётом оптимизации расположения сплайновых узлов// Автометрия. 2013. Т.49. №1. С.26-41.

11. Mandea M., Korte M. Geomagnetic Observations and Models. Springer I AG A Special Sopron Book Series, Vol.5. 1-st Edition. 2011. XV. 343p.

12. Love J. J. Magnetic monitoring of Earth and space// Phvs. Todav.2008. 61. 31-37.

13. Gvishiani A., Lukianova R., Soloviev An. et al. Survey of Geomagnetic Observations Made in the Nothern Sector of Russia and New Methods for Analysing Them// Surveys in Geophysics. 2014. Vol.35. № 5. P.1123-1154. Doi: 10.1107/sl0712-014-9297-8.

14. Gvishiani A., Bogoutdinov S., Soloviev A.A. et al. Automated recognition of spikes in 1 HZ data recordered at the Easter Island magnetic observatory// Earth, Planets and Space. 2012. T. 64. № 9. C. 743-752. DOI: 10.5047/eps.2012.03.004.

15. Gvishiani A.D., Agavan S.M., Bogoutdinov Sh. et al. Recognition of disturbances with specified morphology in time series. Spikes on 1-s magnetograms // Izvestiva. Physics of the Solid Earth. 2012. T. 48. № 5. C. 395-409. DOI: 10.1134/S106935131204009X*

16. Gvishiani A.D., Kihn E., Soloviev A.A. et al. Detection of hardware failures at INTERMASGNET observatories: application of artificial intelligence techniques to geomagnetic records study// Russian Journal of Earth Sciences. 2009. T. 11. № 2. C. lb-6. doi: 10.2205/2009ES000387.

17. Gvishiani A.D., Agavan S.M., Zlotnicki J. et al. Mathematical methods of geoinformatics/ Fuzzy comparision and recognition of anomalies in time series //Cybernetics and Systems Analysis. 2008. T. 44. № 3. C. 309-323. DOI: 10.1007/sl0559-008-9009-9.

18. International Real-Time Magnetic Observatory Network, 2018. http://www. intermagnet.org.

19. Bureau Central de Magnetisme Terrestre,2018. http://www.bcmt.ft\

20. Getmanov V.G., Orlov S.E. A way to Use Local and Spline Models for Estimating the Parametrical Functions of a Nonstationarv Waveform Signals//Pattern Recognition and Image Analysis, 2011. Vol.21, No.4, pp.677-680.

21. Katkovnik V., K.Egiazarian, J.Astola. Local Approximation in Signal and Image Processing. SPIE Publications. 2006. 576p.

22. Seber G.A.F., Lee A.J. Linear Regression Analysis. J.Wilev k, Sons Publications, 2003, 581p. REFERENCES

1. Gvishiani A.D., Luk'vanova R.YU. 2015. "Geoinformatika i nablyudeniya magnitnogo polva Zemli: rossijskij segment" , Fizika Zemli. №2. pp. 3-20. doi:10.7868/S0002333715020040.

2. Getmanov V.G. 2011, Cifrovaya obrabotka nestacionarnyh kolebatel'nyh signalov na osnove lokal'nyh i splajnovyh modelej. M.: Izd-vo NIYAU MIFI. p. 298.

3. Getmanov V.G., Sidorov R.V., Dabagvan R.A. 2015, "Metod fil'tracii signalov s ispol'zovaniem lokal'nyh modelej i funkcij vzveshennogo usredneniva" , Izmeritel'naya tekhnika. №9. pp. 52-57.

4. Getmanov V.G. 1992, "O chastotnom podpoiske v zadache ocenivaniva parametrov kusochno-sinusoidal'nvh funkcij" , Avtometriya. № 2. pp. 93-98.

5. Getmanov V.G., Dabagvan R.A., Sidorov R.V. 2016, "Issledovanie harakteristik geomagnitnvh pul'sacij metodom lokal'nyh approksimacij', Geomagnetizm i aehronomiya. № 2. vol. 56. pp. 209-216.

6. Getmanov V.G. 2009, "Ob algoritme poiska po chastote v zadache ocenivaniva parametrov modelej poligarmonicheskih signalov" , Avtometriya. №3. pp. 83-89.

7. Getmanov V.G., Borzunov G.I. 2015, "Algoritm parallel'nyh vvchislenij dlva zadachi spektral'no- vremennogo analiza na bazisnvh poligarmonicheskih funkcivah" , Informacionnye tekhnologii. vol.21. №6. pp. 456-463.

8. Getmanov V.G. 2013, "Nelinejnava fil'traciya nablvudenij sistemv vektornogo i skalvarnogo magnitometrov" , Izmeritel'naya tekhnika. №6. pp. 51-55.

9. Getmanov V.G. 2016, "Metodv vvchisleniva approksimacionnvh splajnovvh funkcij dlva zadach cifrovoj fil'tracii" , Teoriya i sistemy upravleniya. №5. pp. 45-53.

10. Getmanov V.G. 2013, "Algoritmv vvchisleniva approksimacionnvh splajnovvh funkcij s uchvotom optimizacii raspolozheniva splajnovvh uzlov" , Avtometriya. vol.49. №1. pp. 26-41.

11. Mandea M., Körte M. 2011, Geomagnetic Observations and Models. Springer IAGA Special Sopron Book Series, Vol.5. 1-st Edition. XV. 343 p.

12. Love J. J. 2008, "Magnetic monitoring of Earth and space" , Phys. Today. 61. pp. 31-37.

13. Gvishiani A., Lukianova R., Soloviev An. et al. 2014, "Survey of Geomagnetic Observations Made in the Nothern Sector of Russia and New Methods for Analysing Them" , Surveys in Geophysics. Vol. 35. № 5. pp. 1123-1154. Doi: 10.1107/sl0712-014-9297-8.

14. Gvishiani A., Bogoutdinov S., Soloviev A.A. et al. 2012, "Automated recognition of spikes in 1 HZ data recordered at the Easter Island magnetic observatory" , Earth, Planets and Space. vol. 64. № 9. pp. 743-752. DOI: 10.5047/eps.2012.03.004.

15. Gvishiani A.D., Agavan S.M., Bogoutdinov Sh. et al. 2012, "Recognition of disturbances with specified morphology in time series. Spikes on 1-s magnetograms" , Izvestiya. Physics of the Solid Earth, vol. 48.5. pp. 395-409. DOI: 10.1134/S106935131204009X.

16. Gvishiani A.D., Kihn E., Soloviev A.A. et al. 2009, "Detection of hardware failures at INTERMASGNET observatories: application of artificial intelligence techniques to geomagnetic records study" , Russian Journal of Earth Sciences, vol. 11. № 2. pp. lb-6. doi: 10.2205/2009ES000387.

17. Gvishiani A.D., Agavan S.M., Zlotnicki J. et al. 2008 "Mathematical methods of geoinformatics/ Fuzzy comparision and recognition of anomalies in time series" , Cybernetics and Systems Analysis, vol. 44. № 3. pp. 309-323. DOI: 10.1007/sl0559-008-9009-9.

18. International Real-Time Magnetic Observatory Network, 2018. http://www. intermagnet.org.

19. Bureau Central de Magnetisme Terrestre, 2018. http://www.bcmt.ft\

20. Getmanov V.G., Orlov S.E. 2011 "A way to Use Local and Spline Models for Estimating the Parametrical Functions of a Nonstationarv Waveform Signals" , Pattern Recognition and Image Analysis, Vol.21, No.4, pp.677-680.

21. Katkovnik V., K.Egiazarian, J.Astola. 2006 Local Approximation in Signal and Image Processing. SPIE Publications. 576p.

22. Seber G.A.F., Lee A.J. 2003, Linear Regression Analysis. J.Wilev k, Sons Publications, 581p.

Получено 27.07.2018 Принято в печать 22.10.2018