Введение. В качестве программного инструментария система MATLAB используется в процессе решения задач в самых разных предметных областях [1, 2]. В данной работе рассмотрена задача, результаты решения которой относятся к области геомагнетизма – выделение геомагнитных пульсаций (ГП) из данных, полученных магнитной обсерваторией Lycksele (Швеция, Геологическая служба Швеции, международный IAGA-код LYC), входящей в международную сеть ИНТЕРМАГНЕТ [3–5]. Выбор этой обсерватории обусловлен передачей от нее в ИНТЕРМАГНЕТ достаточно качественных данных измерений состояний геомагнитного поля.

Физически ГП представляют собой изменения напряженности геомагнитного поля в разном диапазоне частот в зависимости от времени и различных факторов. Формально ГП обычно представляют в виде магнитогидродинамических волн в околоземной плазме [6]. Возбуждение таких волн происходит в силу разнообразных физических процессов, которые возникают, например, в ионосфере, магнитосфере или в солнечном ветре. Отсюда следует, что ГП имеют высокую информативность, так как их анализ позволяет получать много информации, например, о параметрах среды в областях генерации ГП, об особенностях развития геомагнитных бурь и суббурь, о потоках высыпающихся частиц и путях их распространения к земной поверхности. Фиксация ГП на поверхности Земли происходит путем регистрации ультранизкочастотных электромагнитных колебаний [6]. Таким образом, можно сделать вывод, что анализ ГП – актуальная задача, связанная с решением как прикладных проблем, например, радиосвязь, диагно- стика текущего состояния околоземной среды и прогноз космической погоды, так и фундаментальных, например, дальнейшее развитие теории ГП [7].

В статье сделан акцент на выделение ГП типа Pc3, то есть ГП с периодом колебаний, попадающим в диапазон 10–45 секунд, и амплитудой магнитной индукции 0,4–20 нТл. Такие ГП относятся к одним из самых распространенных, регистрируемых на земной поверхности. В большинстве случаев ГП Рс3 – это результат резонансных колебаний силовых линий геомагнитного поля в магнитосфере Земли. В работе предложен подход, позволяющий производить выделение ГП Pc3 из данных, полученных в геомагнитных обсерваториях.

Представленные в работе результаты относятся также к области проектирования цифровых фильтров (ЦФ), параметры которых задаются исходя из свойств ГП Pc3. Аналогичные исследования как в России, так и за рубежом начали проводиться достаточно давно. Например, в [8] описаны основы ЦФ, которые использовали тогда в процессе обработки данных измерений с сейсмологических, гравиметрических и магнитометрических датчиков. В настоящее время методы проектирования ЦФ базируются на современных информационных технологиях [9–11]. В последнее время при обработке данных делается акцент на использовании методов технологии Big Data [12]. Необходимо отметить, что публикаций о применении технологий Big Data, реализованных в MATLAB с представлением программного кода, для решения задачи выделения конкретного типа ГП Pc3 в доступных источниках не обнаружено. Вычислительные эксперименты были проведены в системе MATLAB R2023b под управлением Windows 11.

Постановка задачи

Проектирование оптимального рекурсивного полосового ЦФ с бесконечно-импульсной характеристикой (БИХ) необходимо провести при помощи метода Золотарева–Кауэра [13, 14]. Под оптимальным БИХ-фильтром понимается ЦФ минимального порядка при заданных требованиях к амплитудно-частотной характеристике (АЧХ), к полосе пропускания и полосе задержки [13]. Выбор данного ЦФ обусловлен его основными особенностями.

Во-первых, поведение изменений АЧХ как в полосе пропускания, так и в полосе задержки уравновешенное. Амплитуды изменений в любой из этих полос не зависят друг от друга.

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

Во-вторых, порядок ЦФ Золотарева и его сложность минимальные.

Такие ЦФ также относят к эллиптическим. Это связано с тем, что при их описании используют математический аппарат на базе эллиптических функций [13, 15].

Можно отметить, что эти ЦФ называют так из-за дробей Золотарева, применяемых при описании передаточных функций ЦФ. Причем нули передаточных функций расположены так, что минимумы АЧХ находятся в полосе задержки.

Исходя из свойств ГП Pc3 зададим требования к полосе пропускания и к АЧХ в полосе пропускания и полосе задержки проектируемого ЦФ:

f-χ = 1/45 (c) ≈ 20 мГц – левая граничная частота полосы пропускания;

fχ =1/10 (c) ≈ 100 мГц – правая граничная частота полосы пропускания;

amax = 0.001 (дБ) – максимальное допустимое затухание АЧХ в полосе пропускания;

a1min = a2min = 150 (дБ) – максимальное допустимое затухание АЧХ в левой (a1min) и правой (a2min) полосах задержки.

Заключительным требованием при постановке задачи является условие применения реализованных в MATLAB технологий Big Data. Это требование вытекает из большого объема исходных данных: непрерывная съемка данных о состоянии геомагнитного поля по четырем каналам в течение целого года с частотой дискретизации fд = 1 (с). В силу этого необходимо применение специальных средств, таких как Datastore.

Проектирование ЦФ

Зададим параметры для полосового ЦФ, соответствующие принятым условиям. Код на MATLAB можно написать, например, так:

Fs_45_s=1/45; Fs_10_s=1/10; % Диапазон частот Pc3

Fs = 1; % Частота дискретизации 1 секунда

Fpass1 = Fs_45_s; % Левая граница частоты полосы пропускания

Fpass2 = Fs_10_s; % Правая граница частоты полосы пропускания

Fstop1 = 0.99*Fpass1; % Левая граничная частота полосы задержки

Fstop2 = 1.01*Fpass2; % Правая граничная частота полосы задержки

Astop1 = 150; % Величина затухания в левой полосе задержки (дБ)

Astop2 = 150; % Величина затухания в правой полосе задержки (дБ)

Apass = 0.001; % Величина затухания в полосе пропускания (дБ)

Создадим полосовой ЦФ с БИХ-фильтром при помощи следующей команды:

digital_filter_Pc4= designfilt('bandpassiir', ...;

% Полосовой БИХ-фильтр

'StopbandFrequency1', Fstop1, ...; % Левая граничная частота полосы задержки

'PassbandFrequency1', Fpass1, ... ; % Левая граница частоты полосы пропускания

'PassbandFrequency2', Fpass2, ...; % Правая граница частоты полосы пропускания

'StopbandFrequency2', Fstop2, ...; % Правая граничная частота полосы задержки

'StopbandAttenuation1', Astop1, ...; % Величина затухания в левой полосе задержки (дБ)

'PassbandRipple', Apass, ...; % Величина затухания в полосе пропускания (дБ)

'StopbandAttenuation2', Astop2, ...; % Величина затухания в правой полосе задержки (дБ)

'DesignMethod', 'ellip', ...; % Метод Золотарева (эллиптический)

'MatchExactly' ,'passband', ...; % Полосовой ЦФ

'SampleRate', Fs); % Частота дискретизации

В результате выполнения этой команды будет спроектирован оптимальный устойчивый БИХ-фильтр 50-го порядка, имеющий заданные характеристики. На рисунке 1 представлена АЧХ этого ЦФ.

Понятно, что ЦФ можно спроектировать иначе [16]. В этом случае появится дополнительная информация, например, карта нулей  и полюсов передаточных функций у спроектированного ЦФ.

Подготовка данных для проверки  спроектированного цифрового  БИХ-фильтра

Для проверки полученного ЦФ были использованы данные измерений состояний геомагнитного поля в реальном масштабе времени с частотой дискретизации fд = 1 (с), полученные геомагнитной обсерваторией Lycksele за пе- риод с 1 января по 31 декабря 2023 г. Файлы с исходными данными содержат четыре столбца – LYCX, LYCY, LYCZ и LYCF, где LYCX, LYCY, LYCZ – компоненты вектора геомагнитного поля, а LYCF – модуль этого вектора. Все значения приводятся в нанотеслах. Код на MATLAB для доступа к этим данным и их визуализации приведен в [17]. При небольшой модификации, связанной с конкретными вычислительными ресурсами, после выполнения кода появится график исходных данных (рис. 2). Применяя методы визуализации, входящие в Big Data, можно решать множество различных задач, например, максимальное погружение  в данные и (или) обнаружение отклонений и аномалий. На рисунке 2 видно, что в исходных данных имеются выбросы, появившиеся, вероятнее всего, в результате прекращения работы магнитометрических сенсоров.

В таких случаях общепринятой практикой считается занесение числа 99 999 в соответствующую БД. Для дальнейшей работы эти данные необходимо удалить следующим образом:

toDelete_XYZF = (LYC.LYCX == 99999) | (LYC.LYCY == 99999) …

 | (LYC.LYCZ == 99999) | (LYC.LYCF = = 99999);

LYC(toDelete_XYZF,:) = [];

clearvars toDelete_XYZF;

stackedplot(LYC).

Обработанные при помощи такого кода очищенные данные представлены на рисунке 3.

Выделения ГП типа Pc3

На первом шаге необходимо удалить тренд в очищенных данных, например, так: LYC_detrend = detrend(LYC,3); здесь цифра 3 задает тренд в виде кубического сплайна.

Далее пропускаем покомпонентно очищенные данные без тренда через спроектированный ЦФ, учитывая при этом фазовую задержку при помощи функции filtfilt:

LYC.LYCX = filtfilt(digital_filter_Pc3, LYC_detrend.LYCX);

LYC.LYCY = filtfilt(digital_filter_Pc3, LYC_detrend.LYCY);

LYC.LYCZ = filtfilt(digital_filter_Pc3, LYC_detrend.LYCZ);

LYC.LYCF = filtfilt(digital_filter_Pc3, LYC_detrend.LYCF).

В результате выполнения представленного кода в timetable LYC будет содержаться покомпонентно информация о ГП Pc3.

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

Анализ полученных результатов

При помощи методов визуализации можно определить даты, когда происходили магнитные бури. Например, анализируя x-компоненту ГП (LYCX), можно сказать, что 23 и 24 апреля 2023 г. происходила сильная магнитная буря (http://www.swsys.ru/uploaded/image/2024-3/9.jpg).

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

Выводы

Предлагаемый подход дает возможность достаточно эффективно решать задачи предварительной обработки большого объема исходных данных и выделения из них ГП Pc3. Он базируется на проектировании ЦФ и использовании инструментария визуального анализа, позволяющего представлять в графической форме сразу весь набор исследуемых данных. Следует отметить, что в настоящее время при фильтрации геофизических данных также применяют методы спектрального анализа [18, 19].

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

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

1.   Коробейников А.Г. Применение искусственных нейронных сетей в системах автоматического управления магнитной левитацией // Программные продукты и системы. 2022. Т. 35. № 3. С. 452–457. doi: 10.15827/0236-235X. 139.452-457.

2.   Егоров Е.Е. Моделирование работы манипуляционного робота в программном пакете Matlab Robotics Toolbox // Политехнический молодежный журнал. 2020. № 01. С. 1–12.

3.   Коробейников А.Г. Обработка и анализ данных с российского сегмента мировой сети магнитных обсерваторий ИНТЕРМАГНЕТ // Международный журнал гуманитарных и естественных наук. 2018. № 8. С. 91–98.

4.   Вишняков Д.Д., Арутюнян Д.А., Шклярук А.Д. Оценка планетарного индекса геомагнитной активности  по данным стационарных магнитовариационных обсерваторий сети INTERMAGNET с применением собственного программного обеспечения // Гелиогеофизические исследования. 2023. № 38. С. 40–45.

5.   Воробьев А.В., Воробьева Г.Р. Подход к оценке относительной информационной эффективности магнитных обсерваторий сети INTERMAGNET // Геомагнетизм и аэрономия. 2018. Т. 58. № 5. С. 648–652. doi: 10.1134/ s0016794018050164.

6.   Гульельми А.В., Потапов А. С. Частотно-модулированные ультранизкочастотные волны в околоземном космическом пространстве // УФН. 2021. Т. 191. № 5. С. 475–491. doi: 10.3367/UFNr.2020.06.038777.

7.   Mishin V.V., Tsegmed B., Klibanova Yu.Yu., Kurikalova M.A. Burst geomagnetic pulsations as indicators of substorm expansion onsets during storms. JGR Space Phys., 2020, vol. 125, art. e2020JA028521. doi: 10.1029/2020JA028521.

8.   Кулханек О. Введение в цифровую фильтрацию в геофизике; [пер. с англ.]. М.: Недра, 1981. 198 с.

9.   Potthast R., Walter A., Rhodin A. A localized adaptive particle filter within an operational NWP framework. Monthly Weather Review, 2019, vol. 147, pp. 345–362. doi: 10.1175/MWR-D-18-0028.1.

10. van Leeuwen P.J., Kunsch H.R., Nerger L., Potthast R., Reich S. Particle filters for high-dimensional geoscience applications: A review. Quarterly J. of the Royal Meteorological Society, 2019, vol. 145, no. 723, pp. 2335–2365.  doi: 10.1002/qj.3551.

11. Климова Е.Г. Локальный ансамблевый алгоритм усвоения данных для нелинейных геофизических моделей // СибЖВМ. 2023. Т. 26. № 1. С. 27–42. doi: 10.15372/SJNM20230103.

12. Мошак Н.М., Рудинская С.Р., Груздев А.А. Третья платформа информации и Big Data // Наукоемкие технологии в космических исследованиях Земли. 2023. Т. 15. № 4. С. 47–59. doi: 10.36724/2409-5419-2023-15-4-47-59.

13. Солонина А.И. Цифровая обработка сигналов в зеркале MATLAB. СПб: БХВ-Петербург, 2018. 560 с.

14. Будунова К.А., Кравченко В.Ф. Математические методы синтеза частотно-избирательных фильтров // Физические основы приборостроения. 2022. Т. 11. № 1. С. 2–21. doi: 10.25210/jfop-2201-002021.

15. Юлина А.О., Синкевич Г.И. История развития теории эллиптических функций в работах Абеля, Якоби, Вейерштрасса, Сомова // ТВИМ. 2021. № 3. С. 79–92.

16. Коробейников А.Г. Применение MATLAB при проектировании цифровых фильтров для выделения геомагнитных пульсаций типа Pc5 // Программные продукты и системы. 2023. Т. 36. № 2. С. 257–262. doi: 10.15827/ 0236-235X.142.257-262.

17. Коробейников А.Г. Применение методов Big Data для сравнения данных геомагнитных обсерваторий сети INTERMAGNET // Изв. вузов. Приборостроение. 2023. Т. 66. № 12. С. 993–1001. doi: 10.17586/0021-3454-2023-66-12-993-1001.

18. Белашев Б.З. Спектральный анализ данных геофизического мониторинга // Тр. КарНЦ РАН. 2023. № 4.  С. 5–15. doi: 10.17076/mat1768.

19. Лыгин И.В., Брагина А.А., Вишняков Д.Д. Подходы прогноза геомагнитных вариаций // Гелиогеофизические исследования. 2023. № 41. С. 42–53.

References

1. Korobeynikov, A.G. (2022) ‘Applying artificial neural networks in automatic control systems for magnetic levita-tion’, Software & Systems, 35(3), pp. 452–457 (in Russ.). doi: 10.15827/0236-235X.139.452-457.

2. Egorov, E. (2020) ‘Simulating the manipulation robot operation in the Matlab Robotics Toolbox software’, Politechnical Student J., (01), pp. 1–12 (in Russ.).

3. Korobeynikov, A.G. (2018) ‘Processing and analysis of data from the Russian segment of the world network of magnetic observatory INTERMAGNET’, Int. J. of Humanities and Natural Sciences, (8), pp. 91–98 (in Russ.).

4. Vishnyakov, D.D., Arutyunyan, D.A., Shklyaruk, A.D. (2023) ‘Assessment of the planetary index in geomagnetic activity from data of stationary magnetic observatories of the INTERMAGNET network using proprietary software’, Heliogeophysical Research, (38), pp. 40–45 (in Russ.).

5. Vorobev, A.V., Vorobeva, G.R. (2018) ‘Approach to assessment of the relative informational efficiency of INTERMAGNET magnetic observatories’, Geomagnetism and Aeronomy, 58(5), pp. 625–628 (in Russ.). doi: 10.1134/ s0016794018050164.

6. Guglielmi, A.V., Potapov, A.S. (2021) ‘Frequency-modulated ULF waves in near-Earth space’, UFN, 191(5), pp. 475–491 (in Russ.). doi: 10.3367/UFNr.2020.06.038777.

7. Mishin, V.V., Tsegmed, B., Klibanova, Yu.Yu., Kurikalova, M.A. (2020) ‘Burst geomagnetic pulsations as indica-tors of substorm expansion onsets during storms’, JGR Space Phys., 125, art. e2020JA028521. doi: 10.1029/2020JA028521.

8. Kulhánek, O. (1976) Introduction to Digital Filtering in Geophysics. Elsevier Publ., 168 p. (Russ. ed.: (1981) Moscow, 198 p.).

9. Potthast, R., Walter, A., Rhodin, A. (2019) ‘A localized adaptive particle filter within an operational NWP frame-work’, Monthly Weather Review, 147, pp. 345–362. doi: 10.1175/MWR-D-18-0028.1.

10. van Leeuwen, P.J., Kunsch, H.R., Nerger, L., Potthast, R., Reich, S. (2019) ‘Particle filters for high-dimensional geoscience applications: A review’, Quarterly J. of the Royal Meteorological Society, 145(723), pp. 2335–2365. doi: 10.1002/qj.3551.

11. Klimova, E.G. (2023) ‘A local ensemble data assimilation algorithm for nonlinear geophysical models’, SJNM, 26(1), pp. 27–42 (in Russ.). doi: 10.15372/SJNM20230103.

12. Moshak, N.N., Rudinskaya, S.R., Gruzdev, A.A. (2023) ‘Third platform for informatization and big data’, H&ES Research, 15(4), pp. 47–59 (in Russ.). doi: 10.36724/2409-5419-2023-15-4-47-59.

13. Solonina, A.I. (2018) Digital Signal Processing in the MATLAB Mirror. St. Petersburg, 560 p. (in Russ.).

14. Budunova, K.A., Kravchenko, V.F. (2022) ‘Mathematical methods of frequency selective filters synthesis’, Phys. Bases of Instrumentation, 11(1), pp. 2–21 (in Russ.). doi: 10.25210/jfop-2201-002021.

15. Yulina, A.O., Sinkevich, G.I. (2021) ‘The history of the development of the theory of elliptic functions in the works of Abel, Jacobi, Weierstrass, Somov’, TJCSTM, (3), pp. 79–92 (in Russ.).

16. Korobeynikov, A.G. (2023) ‘Applying MATLAB in the design of digital filters for selecting Pc5 geomagnetic pulsations’, Software & Systems, 36(2), pp. 257–262 (in Russ.). doi: 10.15827/0236-235X.142.257-262.

17. Korobeynikov, A.G. (2023) ‘Application of big data methods for comparing data of geomagnetic observatories in the INTERMAGNET network’, J. of Instrument Engineering, 66(12), pp. 993–1001 (in Russ.). doi: 10.17586/0021-3454-2023-66-12-993-1001.

18. Belashev, B.Z. (2023) ‘Spectral analysis of geophysical monitoring data’, Proc. of the Karelian Research Centre of the RAS, (4), pp. 5–15 (in Russ.). doi: 10.17076/mat1768.

19. Lygin, I.V., Bragina, A.A., Vishnyakov, D.D. (2023) ‘Methods for predicting geomagnetic variations’, Heliogeo-physical Research, (41), pp. 42–53 (in Russ.).