Experiment on long-term forecasting of geomagnetic activity based on nonlocal correlations

Full Text

1. Введение

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

Согласно квантовой абсорбционной электродинамике в запутанных состояниях нелокальные корреляции имеют как запаздывающие, так и опережающие компоненты, а их комбинация может вести к кажущимся мгновенным корреляциям [Cramer, 1986; Hoyle and Narlikar, 1995]. В то время как запутывание подсистем, содержащих небольшое число частиц, в настоящее время хорошо изучено, вопрос о макроскопическом запутывании значительно более сложен, его общее теоретическое описание только формируется (хотя такое описание имеется для некоторых статических случаев, например [Amio et al., 2008; Reid et al., 2012; Lee et al., 2015]). В большой серии экспериментов по изучению нелокальных корреляций крупномасштабных случайных диссипативных процессов [Коротаев и др., 2000, 2007, 2015, 2019, 2021; Коротаев и Морозов, 2018; Korotaev, 2011; Korotaev et al., 2005, 2018a,b, 2020, 2021] мы опирались на уравнение макроскопической запутанности, полученное на основе абсорбционной электродинамики Уилера — Фейнмана и ее квантового обобщения. Оно тестировано теоретически (на модели спинового газа [Calsamiglia et al., 2005]) и экспериментально (с магнитосферным источником геомагнитных Dst-вариаций). В соответствии с этим уравнением производства энтропий в случайных диссипативных процессах (пробном процессе в детекторе и процессе-источнике) нелокально коррелируют с симметричным запаздыванием и опережением, которые при диффузионном обмене запутыванием в промежуточной среде могут быть очень велики. Благодаря асимметрии эффективности поглощения запаздывающей и опережающей компонент поля Уилера — Фейнмана промежуточной средой между удаленными процессами [Hoyle and Narlikar, 1995] опережающая корреляция превышает запаздывающую. Асимметрия препятствует их интерференции, приводящей к квазимгновенной корреляции, наблюдаемой в обычных экспериментах с запутанными частицами. Теория детекторов позволяет связать производство энтропии в них с измеряемым сигналом, так же как и в процессах-источниках — с подходящим индексом их активности (хотя на практике проще связывать непосредственно сигнал детектора с этим индексом) [Коротаев и Морозов, 2018; Korotaev, 2011].

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

Основной трудностью в реализации экспериментов является необходимость исключения всех возможных классических локальных воздействий (температура, давление, электромагнитное поле и т.п.) Эта трудность усугубляется тем, что макроскопические нелокальные корреляции проявляются на крайне низких частотах, что требует длительных измерений (год и более). Защита детекторов от крайне низкочастотных помех в условиях обычной лаборатории особенно сложна. Тем не менее эта задача удовлетворительно решена в экспериментальной обсерваторской установке ЦГЭМИ-ИЗМИРАН (Троицк) [Коротаев и Морозов, 2018; Korotaev, 2011]. Что касается природных процессов-источников, то вследствие свойства моногамности квантовых корреляций их число, вызывающее заметную реакцию детектора, не может быть велико, разделение их вкладов в сигнал не вызывает затруднений. Доминирующими оказались глобальные источники — солнечная и геомагнитная активность, вклад регионального источника — синоптической активности — оказался значительно слабее. Подчеркнем, что речь идет именно о нелокальной реакции детекторов на диссипативные процессы в источниках: так, опыты показали отсутствие их отклика в пределах чувствительности на искусственное магнитное поле до 1 мТл в диапазоне 0—1 Гц, а корреляция с потоком космических лучей оказалась значительно слабее, чем непосредственно с индексами геомагнитной и солнечной активности.

Уже в первых сериях измерений длительностью 1.0—2.5 года было подтверждено теоретическое предсказание превышения опережающих нелокальных корреляций над запаздывающими (запаздывающий максимум корреляции мог даже оказаться практически подавленным). Величина опережения глобального максимума оказалась прямо связана с масштабом источника, для глобальной геомагнитной активности, оцениваемой Dst-индексом, она составила десятки-сотни суток, а уровень корреляции оказался достаточным для реализации серий прогностических расчетов Dst с такой заблаговременностью и точностью до 2 нТл [Коротаев и др., 2007; Коротаев и Морозов, 2018; Korotaev, 2011].

Начиная с 2012 г. акцент в этих экспериментах был перенесен на созданную глубоководную установку в озере Байкал. Мощная и спокойная водная толща Байкала предоставляет естественную защиту почти от всех локальных воздействий на электродные детекторы, а возможность их размещения на различных глубинах предоставляет дополнительные возможности для исследований. Единственным недостатком автономной глубоководной установки была необходимость ее подъема на лед в марте каждого года, так что из-за перерыва в работе годовые ряды измерений не могли быть корректно объединены в единый длинный ряд. Установка ЦГЭМИ-ИЗМИРАН в Троицке при этом играла роль опорной обсерватории для разделения вкладов в сигналы детекторов глобальных и региональных процессов. Благодаря этому в сигналах детекторов байкальской установки удалось выделить опережающую реакцию на сильные землетрясения и на макротурбулентность в деятельном слое, причем для последних выполнить серии успешных прогностических расчетов с заблаговременностью до полутора месяцев. Опережающие нелокальные корреляции с солнечно-геомагнитной активностью в сигналах байкальских детекторов также всегда наблюдались, и возможность их прогноза многократно продемонстрирована простым сдвигом соответственно фильтрованным рядам, но их полноценный количественный прогноз был невозможен. Дело в том, что любые самообучающиеся прогностические алгоритмы требуют длины обучающего интервала существенно большей заблаговременности прогноза, а оптимальная заблаговременность (определяемая положением глобального максимума корреляционной функции) достигала для солнечно-геомагнитной активности 8—10 месяцев.

Целью данной работы являются решение задачи использования накопленного на экспериментальной установке ЦГЭМИ-ИЗМИРАН материала для тестового прогноза глобальной геомагнитной активности и сравнение различных прогностических алгоритмов.

2. Эксперимент

Установка ЦГЭМИ-ИЗМИРАН, детально описанная в работах [Коротаев и Морозов, 2018; Korotaev, 2011] включает электродный и фотокатодный детекторы макроскопических нелокальных корреляций и аппаратуру для сопровождающих измерений. Данные электродного детектора более надежны, и с ними получены более длинные ряды наблюдений, поэтому именно они использованы в настоящей работе.

Этот детектор включает подобранную пару метрологических морских серебряных-хлорсеребряных слабополяризующихся электродов ХД-5.519.00, выпущенных единственной малой серией фирмой “Эталон” (Киев). Эти электроды были разработаны для прецизионных измерений постоянного электрического поля в море и являются лучшими в мире по нечувствительности собственных потенциалов к физико-химическим условиям окружения. Фактически единственным известным влияющим на них фактором является температура (средний температурный коэффициент — 0.04 мВ/К), что может быть учтено ее контролем. Упрощенная схема устройства электродного детектора показана на рис. 1 (сложное внутреннее устройство электродов не показано). Температура контролируется тремя датчиками, два из которых (t₁ и t₂) установлены у каждого электрода, а третий (t_e) — на расстоянии 1 м, характеризуя общую температуру в лаборатории на высоте детектора.

 

Рис. 1. Схема устройства детектора [Korotaev et al., 2019]: Е — электроды (сложное внутреннее устройство не показано), разность потенциалов которых (U) является сигналом детектора; NaCl — электролит (3% водный раствор); V — корпус детектора (эбонит); t₁,t₂.t_e — термодатчики. Расстояние между электродами — 30 см.

 

Влияние температуры было изучено предварительно на основе многомесячных непрерывных измерений естественных вариаций: средний температурный коэффициент — 0.040 мВ/К, запаздывание (обусловленное стабилизирующим действием водного электролита) — 90 мин. Сигнал детектора U и температура измерялись с частотой 10 Гц и осреднялись за каждую минуту.

Указанный температурный коэффициент является именно средним, само это понятие является условным из-за нелинейности и частотной зависимости, в итоге — нестационарности влияния температуры на U при длительных измерениях. Поэтому учет влияния температуры велся для каждого отсчета методом текущей регрессии. Коэффициент регрессии вычислялся (с учетом запаздывания) за 28-суточный интервал, предшествующий индивидуально исправляемому отсчету. Итоговая относительная погрешность измерений U составляет не более 0.01%.

В процессе длительных исследований с установкой она неоднократно усовершенствовалась, поэтому в измерениях происходили существенные перерывы. Для обработки выбран самый длинный имеющийся непрерывный ряд наблюдений: с 25.09.2009 г. по 15.04.2017 г.

Как отмечено во Введении, опережающие нелокальные корреляции наблюдаются только для чисто случайных процессов. Длины рядов, использованных в предыдущих экспериментах по прогнозу Dst на основе нелокальных корреляций с U [Коротаев и др., 2007; Коротаев и Морозов, 2018; Korotaev, 2011], не позволяли достоверно охватить диапазон квазипериодов более года. Поэтому полосы квазипериодов случайных колебаний Dst между периодами, детерминированными обращением Земли вокруг Солнца и собственным вращением Солнца и их гармониками, были весьма узки. Длина ряда U, используемого в настоящем эксперименте (7.5 лет), позволяет охватить диапазон от 1 до 4 лет, свободный как от внутригодовых периодичностей, так и от гармоник 11-летнего цикла, т.е. диапазон, где случайная компонента в изменчивости Dst преобладает.

3. Обработка данных

Использовались ряды среднесуточных значений U и Dst. Для расчета корреляционной функции и функции независимости ряд Dst взят с началом на год раньше начала ряда U и с концом на год позже U. Это обеспечивает равноточные определения функций при всех положительных и отрицательных сдвигах времени. Сдвиги обоих знаков нужны для того чтобы убедиться, что опережающий корреляционный максимум является глобальным, что исключает предположение о воздействии геомагнитной активности на сигнал детектора каким-либо локальным механизмом.

Данные подвергались предварительной широкополосной фильтрации, границы полосы подбирались эмпирически для достижения максимальной корреляции. Оптимальная фильтрация соответствует полосе периодов 1307 > T > 365 сут.

По фильтрованным данным в диапазоне сдвигов времени −365 < τ < +365 сут рассчитывалась корреляционная функция U и Dst и функция независимости U от Dst. Последняя определяется через условную и безусловную шенноновские энтропии переменных: $i_{U|Dst}=S\left(U\right|Dst)/S\left(U\right)$ , она ведет себя, грубо говоря, обратно модулю корреляции, но, в отличие от нее, адекватна при любой нелинейности связи переменных. Результаты расчета показаны на рис. 2. Глобальный максимум корреляции $r=0.97\pm 0.01$  при опережении U относительно Dst на 329 сут. Ему соответствует глобальный минимум независимости $i_{U|Dst}={0.37}_{-0.01}^{+0.02}$.

 

Рис. 2. Корреляционная функция r U и Dst и функция независимости $i_{U|Dst}$ τ — сдвиг времени Dst относительно U, сут.

 

4. Прогностические расчеты и результаты

Простейшим способом демонстрации возможности прогноза Dst по U является сдвиг реализаций на время, соответствующее опережению глобального максимума корреляции (рис. 3). Видно, что кривая U действительно в сглаженном виде предсказывает кривую Dst с такой большой заблаговременностью. Но это именно демонстрация возможности прогноза, но не настоящий прогноз, потому что здесь нет количественного предсказания величин Dst — на рисунке кривые просто масштабированы так, чтобы их подобие было наглядным, количественная оценка связи, конечно, возможна, но лишь post factum. Кроме того, кажущаяся гладкость кривой U является лишь следствием масштаба, т.е. связь U с Dst не определяется единственным и постоянным коэффициентом. Для настоящего прогноза требуется сначала нахождение зависимости Dst от U (в рамках тех или иных допущений), а затем ее использование для прогноза.

 

Рис. 3. Демонстрация возможности прогноза по U путем сдвига фильтрованных рядов на величину опережения глобального максимума корреляции (329 сут). Верхняя ось времени соответствует U, нижняя — Dst.

 

Самым простым самообучающимся прогностическим алгоритмом является алгоритм текущей регрессии. На некотором первоначальном обучающем интервале сдвинутых на величину оптимальной заблаговременности (в данном случае 329 сут) рядов рассчитываются параметры линейной регрессии. Далее ряд U сдвигается вперед, и прогнозируется единственное значение Dst с той же заблаговременностью. Затем обучающий интервал подтягивается на один дискрет (в данном случае 1 сут), прогнозируется следующее значение Dst и т.д. В итоге получается прогностическая серия Dst c фиксированной заблаговременностью. Этот алгоритм учитывает нестационарность связи в рамках марковского и кусочно-линейного приближения. Результат его применения показан на рис. 4. При расчете оптимальная длина обучающего интервала составила 500 сут, поэтому за счет выброса первичного интервала прогностическая серия на этом рисунке начинается на 500 сут позже, чем на рис. 3. Прогностическая кривая на рис. 4 похожа на кривую U на рис. 3, она также сглажена по сравнению с фактической, хотя в целом ближе к ней. Но главное, что здесь представлен количественный прогноз, а его точность достаточна для всех практических целей.

 

Рис. 4. Прогноз Dst методом текущей регрессии с фиксированной заблаговременностью 329 сут (тонкая линия) в сопоставлении с фактической кривой (жирная линия). Среднеквадратичная погрешность прогноза — 0.99 нТл.

 

Алгоритм текущей импульсной характеристики (разработанный и использованный еще в первом успешном опыте прогноза солнечной и геомагнитной активности на основе нелокальных корреляций [Коротаев и др., 2007]) также основан на кусочно-линейном приближении, но учитывает немарковский характер процесса. На обучающем интервале $[t_1,t_2]$ вычисляется импульсная переходная характеристика $g\left(\tau \right),$, которая связывает сигнал детектора U и прогнозируемый индекс Dst путем решения уравнения:

$Dst\left(t\right)=\underset{t_1}{\overset{t_n}{\int }}g(\tau \left)U\right(t-\tau )d\tau .$ (1)

Решение интегрального уравнения (1) (обратная задача) в дискретной форме сводится к решению системы линейных уравнений $\{\mathbf{Dst}=\mathbf{UK}\}$. Компоненты вектора K эквивалентны коэффициентам множественной регрессии (в случае собственного распределения). Число уравнений $n$  равно заблаговременности прогноза. U— квадратная матрица $n\times n$, строки которой формируются из значений сигнала детектора на обучающем интервале. Первая строка содержит значения с временным индексом от 1 до $n$, вторая — от 2 до n n + 1 и т.д. Последовательные значения Dst ставятся в соответствие каждой строке матрицы. Система решается методом Гаусса. Устойчивость результатов достигается путем оптимальной регуляризации. Практически опережение выбирается равным ожидаемому среднему положению главного корреляционного максимума. Полный обучающий интервал Dst заканчивается последним наблюденным значением, а U — предшествующим на $\Delta t$.

Найденная таким путем импульсная переходная характеристика используется на следующем интервале данных сигнала детектора, заканчивающемся последним измеренным значением, и прогнозируется одно значение Dst с заблаговременностью $\Delta t$ . Таким образом, решается прямая задача (1) по интервалу U, заканчивающемуся последним наблюденным значением. На следующем шаге обучающий интервал подтягивается вперед на один дискрет, процедура повторяется, и прогнозируется следующее значение. Это минимизирует влияние нестационарности и нелинейности.

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

На рис. 5 представлен результат прогноза этим методом. Длина обучающего интервала составила 658 сут (минимально допустимая методом длина в данном приложении), поэтому тестовый ряд получился короче, чем на рис. 4, но детальность прогноза значительно лучше. Среднеквадратичная погрешность 0.40 нТл. Таким образом, учет немарковского характера процесса увеличил точность прогноза в 2.5 раза.

 

Рис. 5. Прогноз Dst методом текущей импульсной переходной характеристики с фиксированной заблаговременностью 329 сут (тонкая линия) в сопоставлении с фактической кривой (жирная линия). Среднеквадратичная погрешность прогноза — 0.40 нТл.

 

Для прогноза Dst по U была также использована нейронная сеть, что позволило отказаться от кусочно-линейного приближения связи переменных. Нейросетевой метод обратного распространения ошибки с успехом применяется для прогнозирования в самых разных областях геофизики. Методика прогнозирования сейсмических событий по низкочастотным сигналам была представлена в работах [Попова и др., 2016; Popova et al., 2018]. Прогноз магнитуды землетрясения за два дня до его возникновения по изменениям параметров геомагнитного поля, атмосферы и температуры поверхности Земли был выполнен в работе [Suratgar et al., 2008]. Предсказание магнитуды сейсмического события по сейсмическим данным, электрическому сигналу и магнитуде предыдущего землетрясения рассматривалось в работе [Dutta, 2011].

Алгоритм обратного распространения ошибки, как правило, использует трехслойную сеть, состоящую из входного слоя нейронов, на который подаются входные сигналы, выходного слоя нейронов и скрытого слоя нейронов между ними. Нейроны каждого слоя связаны с нейронами последующего слоя, и каждой связи приписывается весовой коэффициент $W_{ij}$, который соединяет i- и j- нейроны. Решение задачи на такой структуре реализуется поэтапно. На первом шаге нейронная сеть обучается на известных примерах. Обучение заключается в подборе таких значений весовых коэффициентов для всех слоев сети, чтобы при заданном входном векторе из обучающей выборки получить на выходном слое сети соответствующий ему выходной вектор с требуемой точностью.

Алгоритм обратного распространения определяет стратегию подбора весов сети с применением градиентного метода оптимизации. Основу этого метода составляет целевая функция, формулируемая в виде квадратичной суммы разностей между фактическими $u_i$ и ожидаемыми $u_i^t$ значениями сигналов на выходном слое сети:

$Er=\sum _p\sum _i{\left(u_i-u_i^t\right)}^2.$ (2)

Суммирование (2) выполняется по всем обучающим примерам $p$ и нейронам $i$ выходного слоя. Значение  формируется на выходном нейроне как результат прохождения сигнала от входного слоя к выходному.

Обучение состоит в вычислении матрицы связей нейронов $W_{ij}$ для каждой пары слоев в целях минимизации ошибки (2). Обучение начинается с небольших случайных значений, которые затем устанавливаются в процессе минимизации ошибки методом наискорейшего спуска:

$\Delta W_{ij}^{\left(n\right)}=-\alpha \frac{\partial Er}{\partial W_{ij}}+\beta \Delta W_{ij}^{(n-1)},$ (3)

где $\Delta W_{ij}^{\left(n\right)}$ — приращение матрицы связей на шаге $n$; $\Delta W_{ij}^{(n-1)}$ — приращение на предыдущем шаге; $\alpha$ и $\beta$ — внутренние параметры нейросети. Процедура (3) выполняется для всего обучающего ряда и заканчивается при достижении некоторого задаваемого пользователем порога точности $eps$, когда $Er$< $eps$.

Обучение нейронной сети для последующего прогнозирования временного ряда Dst-индекса геомагнитной активности выполнялось следующим образом. В качестве входов нейросети выбирались сигналы детектора U из временного окна t(1) – t(n), сдвинутые относительно выхода — примера единственного значения Dst на величину оптимальной заблаговременности. Таким образом формировался первый пример из обучающей выборки. Затем выбранное временное окно для значений U сдвигалось на один шаг вперед, и ему ставилось в соответствие сдвинутое на шаг вперед значение Dst и т.д. Исследования показали, что для выбранного временного окна, размер которого соответствовал 200 сут, длина обучающего интервала для корректного прогноза составила 642 сут. После обучения нейросеть используется для следующего прогноза со сдвигом в несколько дней, продолжая обучаться на уже пройденных временных участках.

Прогностическая кривая представлена на рис. 6. Среднеквадратичная погрешность 0.29 нТл. Таким образом, полный учет нелинейности связи (отказ от кусочно-линейного приближения) ведет к увеличению точности прогноза примерно в 1.5 раза.

 

Рис. 6. Прогноз Dst методом текущей нейросети с фиксированной заблаговременностью 329 сут (тонкая линия) в сопоставлении с фактической кривой (жирная линия). Среднеквадратичная погрешность прогноза — 0.29 нТл.

 

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

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

Реализация такого прогноза Dst-индекса продемонстрирована тремя методами на длительной временной серии с фиксированной заблаговременностью около 11 месяцев (329 сут), отвечающей положению глобального максимума корреляционной функции. Наиболее простой метод текущей регрессии, учитывающий нестационарность нелокальной связи источника и пробного процесса в детекторе в рамках кусочно-линейного приближения, продемонстрировал среднеквадратичную погрешность 0.99 нТл. Метод текущей импульсной характеристики, учитывающий в рамках того же приближения немарковский характер процессов, показал погрешность 0.40 нТл. Метод текущей нейросети, допускающий, кроме того, любую нелинейность связи, показал наименьшую погрешность 0.29 нТл. Таким образом, точность прогноза Dst на основе нелокальных корреляций достаточна для практических целей, в том числе с применением наиболее простого алгоритма.

About the authors

S. M. Korotaev

Schmidt Institute of Physics of the Earth, RAS (IPE RAS)

Author for correspondence.
Email: korotaev@gemrc.ru

Geoelectromagnetic Research Centre (GEMRC)

Russian Federation, Moscow, Troitsk

V. O. Serdyuk

Schmidt Institute of Physics of the Earth, RAS (IPE RAS)

Email: korotaev@gemrc.ru

Geoelectromagnetic Research Centre (GEMRC)

Russian Federation, Moscow, Troitsk

I. V. Popova

Schmidt Institute of Physics of the Earth, RAS (IPE RAS)

Email: korotaev@gemrc.ru

Geoelectromagnetic Research Centre (GEMRC)

Russian Federation, Moscow, Troitsk

J. V. Gorohov

Pushkov Institute of Terrestrial Magnetism, Ionosphere and Radio Wave Propagation, RAS IZMIRAN)

Email: korotaev@gemrc.ru
Russian Federation, Moscow, Troitsk

E. O. Kiktenko

Schmidt Institute of Physics of the Earth, RAS (IPE RAS)

Email: korotaev@gemrc.ru

Geoelectromagnetic Research Centre (GEMRC)

Russian Federation, Moscow, Troitsk

D. A. Orekhova

Schmidt Institute of Physics of the Earth, RAS (IPE RAS)

Email: korotaev@gemrc.ru

Geoelectromagnetic Research Centre (GEMRC)

Russian Federation, Moscow, Troitsk

References

Supplementary files