Восстановление векторного магнитного поля по однокомпонентным данным с использованием эволюционного алгоритма

Полный текст

1. ВВЕДЕНИЕ

Измерение векторного магнитного поля необходимо для ряда прикладных задач, например, для нужд навигации и направленного бурения [Kaji et al., 2019; Buchanan et al., 2013; Колесова, 1985].

Полное магнитное поле описывается трехмерным вектором (H_x, H_y, H_z), модуль которого равен $T=\sqrt{H_x^2+H_y^2+H_z^2}$. Однако, в зависимости от типа магнитометра, не всегда удается провести измерения всех трех компонент поля. Например, множество измерений локальных магнитных полей в геофизике проводится с помощью протонного магнитометра. Такой магнитометр почти не требует специальной ориентации в пространстве, но с его помощью измеряется только модуль полного магнитного поля T [Liuet al., 2022]. Однако для приложений в навигации и направленном бурении необходима карта именно векторного магнитного поля [Park and Choi, 2020; Kaji et al., 2019].

Альтернативой проведению прямых измерений векторного магнитного поля являются методы расчета, продолжения и сепарации магнитных полей по данным одной компоненты [Arturi et al., 2004; Zuo et al., 2020; Baniamerian et al., 2020; Munschy and Fleury, 2011; Pilkington and Boulanger, 2017]. Один из таких методов был предложен в работе [Lourenco and Morrison, 1973], где компоненты поля описывались двойными рядами Фурье, коэффициенты которых связаны через потенциал аномального магнитного поля. Однако эта задача рассматривалась в предположении, что магнитные моменты всех источников ориентированы только в направлении главного магнитного поля Земли на всей площади наблюдений.

Для решения аналогичных задач по поиску и определению параметров залегающих магнитных источников, широко используются методы из семейства мягких вычислений [Ibrahim, 2016]. К их числу относятся эволюционные и генетические алгоритмы [Montesinos et al., 2016; Maier et al., 2019; Kaftan, 2017; Ding et al., 2022], алгоритм имитации отжига [Biswas and Acharya, 2016], оптимизация роем частиц [Pace et al., 2021; Essa and Elhussein, 2020], которые уже успешно применялись для инверсии потенциальных полей в задачах разведочной геофизики.

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

Алгоритм был протестирован на синтетическом примере карты магнитного поля, которая создается 50 случайно ориентированными диполями, распределенными случайно в 3D пространстве под поверхностью наблюдений. Далее, на основе данных спутниковой съемки геомагнитного поля, представленных в моделях IGRF, которая использовалась как нормальное поле [Alken, 2021] и EMM [The National Centers for Environmental Information, 2018] в качестве полного поля, вычислялось аномальное поле и с помощью разработанного эволюционного алгоритма были восстановлены H_x, H_y, H_z-компоненты магнитного поля по вычисленному на высоте 100 м над уровнем моря в прямоугольной области площадью 300×300 км² над полуостровом Ямал аномальному полю.

2. ОПИСАНИЕ АЛГОРИТМА

В качестве источников аномального магнитного поля в данной работе используются точечные диполи, так как магнитное поле намагниченного тела на расстояниях, превышающих его размеры, эквивалентно полю точечного диполя. Поле точечного магнитного диполя определяется известной формулой [Яновский, 1978]

$\overrightarrow{H}\left(\overrightarrow{r}\right)=\frac{3\overrightarrow{R}(\overrightarrow{m},\overrightarrow{R})}{{\left|\overrightarrow{R}\right|}^5}-\frac{\overrightarrow{m}}{{\left|\overrightarrow{R}\right|}^3}$, (1)

где $\overrightarrow{R}$ – вектор из точки нахождения диполя в точку измерения, $\overrightarrow{m}=(m_x,m_y,m_z)$ − магнитный момент диполя в декартовой системе координат.

Общая геометрия рассматриваемой задачи показана на рис. 1a. Пусть векторное магнитное поле необходимо восстановить по данным измерения одной его компоненты на плоскости L_x × L_y, на высоте h над земной поверхностью. Восстановление поля проводится с помощью набора диполей S = {d_i, i = 1…N_d}, залегающих на глубине не более L_z от поверхности земли. Каждый точечный диполь d_i из набора S характеризуется шестью параметрами d_i = {x_i, y_i, z_i, M_i, I_i, D_i}, где x_i, y_i, z_i – положение точечного i-го диполя в декартовой системе координат, M_i – модуль магнитного момента i-го диполя, а I_i и D_i – склонение и наклонение магнитного момента i-го диполя, так что компоненты магнитного момента диполя равны m_xi = M_i cos I_i cos D_i, m_zi = M_i sin I_i.

 

Рис. 1. (a)− Схематичное изображение набора пробных диполей и создаваемого ими аномального магнитного поля, иллюстрирующее геометрию задачи; (б) − блок-схема эволюционного алгоритма для расчета аномального магнитного поля с помощью набора точечных дипольных источников.

 

Решение обратной задачи по восстановлению векторного поля происходит путем подбора координат и магнитных моментов диполей из набора S так, чтобы создаваемое ими однокомпонентное магнитное поле H^S(x, y, h) совпадало с принятой точностью с измеренным однокомпонентным полем H^exp(x, y, h). Пусть, например, по условию задачи известна только вертикальная компонента магнитного поля, H_z(x, y, h). Тогда, после того, как оптимальный набор диполей найден, поперечные компоненты поля в рассматриваемой плоскости, H_x(x, y, h) и H_y(x, y, h), вычисляются с использованием формулы (1).

Поиск оптимального распределения магнитных диполей в данной работе реализуется с помощью эволюционного алгоритма [Maier et al., 2019], который имитирует процесс естественного отбора в биологической системе. Потенциальные решения оптимизационной задачи, то есть в данном случае множество наборов диполей S^j, рассматриваются как особи некоторой биологической популяции. К популяции последовательно применяются эволюционные операторы скрещивания и мутации. Затем, качество каждой особи оценивается с помощью специальной fitness-функции. Наконец, к популяции применяется оператор селекции, в результате которого лучшие особи популяции отбираются и переходят на следующую итерацию алгоритма, после чего цикл алгоритма повторяется. Работа алгоритма продолжается до тех пор, пока значение fitness-функции не уменьшится до требуемой величины, и оптимальное решение в смысле данной fitness-функции не будет найдено. Блок схема эволюционного алгоритма показана на рис. 1б. Ниже приведены описания операторов эволюционного алгоритма.

2.1. Инициализация

Исходными данными при работе алгоритма являются известные размеры области сканирования магнитного поля на плоскости, L_x и L_y, высота сканирования h, количество N_st точек, в которых проводились измерения аномального магнитного поля, их координаты и значения известной компоненты магнитного поля, например, компоненты H_z.

На этапе инициализации алгоритма задается максимальная глубина L_z залегания точечных диполей и их скалярные магнитные моменты M. Потенциальными решениями задачи являются наборы диполей S_j, j = 1…N_s. В каждом наборе имеется N_d диполей. В реализованном варианте алгоритма N_s = 10. Параметры диполей d_i = {x_i, y_i, z_i, I_i, D_i}, i = 1…N_d, то есть их положения и ориентация, на этапе инициализации алгоритма задаются случайно, а оптимальное количество диполей N_d и их магнитные моменты M выбираются с помощью специальной процедуры.

2.2. Скрещивание

К случайно созданной популяции сначала применяется оператор скрещивания. Оператор скрещивания случайным образом выбирает из популяции два различных набора S_a и S_b и создает новый набор S_c из случайно выбранных диполей наборов S_a и S_b, при этом число диполей N_d в наборе сохраняется. Получившийся новый набор S_c добавляется к общей популяции. В реализованной версии алгоритма с помощью операции скрещивания в общую популяцию добавляется 5 дополнительных наборов S, так что общее число наборов в популяции становится равным N_s = 15.

2.3 Мутация

Оператор мутации выбирает произвольный набор S_m из популяции S_j, j = 1…N_s, где N_s = 15, и затем вносит малые случайные изменения в параметры диполей этого набора. Оператор мутации также применяется к популяции 5 раз, в результате чего общее число наборов популяции становится равным N_s = 20.

2.4 Селекция

Оператор селекции производит оценку качества отдельных членов популяции с помощью расчета для каждого члена популяции S_j, j = 1…20 целевой функции

$f_j=\frac{1}{N_{st}}\sum _{x,y}\frac{\left|H^{\exp }(x,y,h)-H^{S_j}(x,y,h)\right|}{⟨H^{\exp }(x,y,h)⟩}$. (2)

В формуле (2) H^exp(x, y, h) есть измеренная компонента аномального магнитного поля в точках измерения (x, y, h), $⟨H^{\exp }(x,y,h)⟩=\sum _i^{N_{st}}\left|H^{\exp }\left(x_i,y_i,h\right)\right|/N_{st}$ − средняя компонента поля на всей области измерений, H^S(x,y,h) есть соответствующая компонента поля пробных диполей из набора S в тех же точках, N_st – общее число точек измерений, а суммирование производится по узлам сетки, на которой проводились измерения поля в плоскости над источниками. Целевая функция (2) в худшем случае принимает значения f ≈ 1, но может уменьшаться до f = 0, если истинное и пробное магнитные поля полностью совпали. Проведенные тестовые расчеты показывают, что, как правило, векторное магнитное поле восстанавливается с хорошей точностью при значениях целевой функции f < 0.2.

Наконец, алгоритм выбирает из общей популяции 10 лучших наборов диполей с наименьшими значениями целевой функции, а остальные наборы удаляет. Ко вновь образованной популяции снова последовательно применяются операторы скрещивания, мутации и селекции до тех пор, пока для наилучшего набора диполей fitness-функция не достигнет достаточно малых значений, порядка f₀ = 0.05 – 0.2.

Для достижения указанной точности алгоритму обычно требуется 50 – 500 итераций в простом случае, если восстанавливаемое магнитное поле создается 5 – 10 дипольными источниками и 10³ – 10⁵ итераций, если количество диполей N_d > 10.

2.5 Усреднение решений

С помощью эволюционного алгоритма можно последовательно получить несколько различных оптимальных распределений диполей с близкими и достаточно малыми значениями fitness-функции. Результирующее аномальное магнитное поле в данной работе строится путем суперпозиции, то есть объединения всех найденных лучших решений. При этом модули магнитных моментов всех диполей при расчете компонент магнитного поля уменьшаются в соответствующее число раз, так что $M_{final}^i=M^i/N_{best}$, где N_best есть число найденных оптимальных наборов диполей с близкими значениями fitness-функции. Результаты проведенных расчетов показывают, что операция суперпозиции позволяет заметно уменьшить итоговое значение fitness-функции при расчетах с большим количеством диполей, а также при восстановлении аномальных магнитных полей сложной конфигурации.

3. СЛУЧАЙНО РАСПРЕДЕЛЕННЫЕ ДИПОЛИ

Разработанный алгоритм был протестирован путем расчета поперечных компонент магнитного поля H_x, H_y по известной вертикальной компоненте поля H_z, заданной на некоторой плоскости. Тестовое магнитное поле создается 50 точечными магнитными диполями с фиксированной величиной магнитного момента M = 3∙10^-3 А∙м². Направления же магнитных моментов и координаты диполей при этом задаются случайно в области с размерами L_x = 100 см, L_y = 100 см и L_z = 10 см, высота сканирования магнитного поля h = 10 см.

При восстановлении векторного магнитного поля по заданным значениям H_z компоненты, необходимо определить количество пробных диполей N_d, их средний момент M, а также положения и направления моментов пробных диполей d_i = {x_i, y_i, z_i, I_i, D_i}, i = 1…N_d в указанной области. Количество варьируемых параметров задачи при этом крайне велико, поэтому на первом этапе определяются оптимальные значения величин M и N_d с помощью следующей процедуры.

Сначала рассчитывается среднее значение заданной компоненты магнитного поля на известной плоскости, $⟨H_z⟩=\sum _i^{N_{st}}\left|H_z^{}\left(x_i,y_i,h\right)\right|/N_{st}$. Эту же величину можно оценить, создавая в области поиска случайные распределения пробных диполей с различными значениями M и N_d. Оптимальные значения указанных величин можно получить, минимизируя функцию невязки $f_{er}=\left|⟨H_z^{\exp }⟩-⟨H_z⟩\right|/⟨H_z^{\exp }⟩$. Распределение функции невязки для диапазона значений магнитного момента M = 10^-4 – 10^-3 А∙м² и количества диполей N_d = 1 − 100 для заданного тестового магнитного поля $⟨H_z^{\exp }⟩$ показано на рис. 2.

 

Рис. 2. Распределение функции невязки для различных значений магнитных моментов диполей, M, и числа магнитных диполей, N_d. Римскими цифрами I − III обозначены области параметров M и N_d, для которых проводились расчеты с помощью эволюционного алгоритма.

 

На рис. 2 видно, что оптимальным значениям параметров M и N_d соответствует значение функции невязки f_er < 0.2. Стоит отметить, что полный расчет диаграммы на рис. 2 занимает около минуты на персональном компьютере. Этот прием можно использовать для предварительной оценки оптимальных значений магнитных моментов и количества дипольных источников, необходимых для восстановления конкретного векторного магнитного поля. На основе такой диаграммы можно выбрать несколько лучших пар значений M и N_d, для которых и проводится дальнейший расчет с помощью эволюционного алгоритма. В данной работе для проведения дальнейших расчетов были выбраны области I–III, показанные рис. 2, которым соответствуют различные M и N_d, приведенные в табл. 1.

 

Таблица 1. Значения параметров M и N_d для проведения оптимизационной процедуры, а также финальное значение fitness-функции

	Магнитный момент, M, А∙м2	Число диполей, Nd	Финальное значение fitness-функции
Область I	5.5∙10-3	20	0.045
Область II	3∙10-3	50	0.05
Область III	2.3∙10-3	80	0.05

 

Из табл. 1 видно, что оптимальная область II соответствует заданным значениям M = 3∙10^-3 А∙м² и N_d = 50. Однако, как показывает рис. 2, оптимальными также являются область I, включающая меньшее число диполей, но с большим магнитным моментом, а также область III, включающая большее число диполей с меньшим магнитным моментом. Расчеты показывают, что финальные значения fitness-функции для областей I–III совпадают с хорошей точностью и достигают значений около 0.2. Это показывает, что наборы диполей I – III фактически являются эквивалентными.

Результаты восстановления векторного магнитного поля с параметрами M и N_d из области III показаны на рис. 3. В данном тестовом расчете магнитное поле задавалось в 400 дискретных точках, которые показаны на рис. 3а. Расстояние между соседними точками составляет 1 см. Заданные значения вертикальной H_z компоненты магнитного поля в указанных точках плоскости поступали на вход эволюционного алгоритма и использовались при вычислении fitness-функции по уравнению (2).

 

Рис. 3. (а, б, в) − Заданное распределение компонент H_x, H_y, H_z магнитного поля;

(г, д, е) − найденные с помощью эволюционного алгоритма распределения тех же компонент H_z, H_x, H_y, по данным одной вертикальной компоненты H_z для случая значений M = 2.3∙10^-3 А∙м² и N_d = 80 (область III на рис. 2); (ж, з, и) – распределение невязки для компонент найденного магнитного поля.

 

На рис. 3 каждая масштабная шкала соответствует одному из распределений H_z, H_x или H_y компонент магнитного поля. Как показывают рис. 3а–в, абсолютная величина каждой компоненты магнитного поля достигает значений порядка 350–450 нТ. На вход эволюционного алгоритма поступали только данные о распределении перпендикулярной H_z компоненты магнитного поля, рис. 3а. Поперечные же заданные компоненты поля H_x и H_y приведены на рис. 3б и 3в для сравнения с восстановленными компонентами, показанными на рис. 3д и 3е, соответственно.

Найденные с помощью эволюционного алгоритма контуры магнитных полей H_z, H_x и H_y на рис. 3г и 3е хорошо совпадают с соответствующими известными контурами магнитного поля, представленными на рис. 3а–в.

На рисунках 3ж–и представлены распределения невязки для трёх компонент магнитного поля. Здесь видно, что на основной площади измерений удается достичь локальных значений невязки, много меньших 0.05. Однако из-за влияния локальных краевых эффектов, рис. 3з–и, общая невязка увеличивается, f_er = 0.05. Если же исключить края из расчета ошибки, то значение fitness-функции в данном примере f_er = 0.015.

Подчеркнем, что результаты, представленные на рис. 3г и рис. 3е, были получены для параметров M и N_d из области III. Однако близкие распределения полей были получены также для расчетов с параметрами M и N_d из областей I и II на рис. 2.

4. ВЫЧИСЛЕНИЕ ВЕКТОРНОГО ПОЛЯ ПО ИЗМЕРЕННЫМ ДАННЫМ АНОМАЛЬНОЙ СОСТАВЛЯЮЩЕЙ ГЕОМАГНИТНОГО ПОЛЯ

С помощью разработанного эволюционного алгоритма можно также рассчитать векторное магнитное поле по данным аномальной составляющей полного геомагнитного поля. Для получения модульных данных аномального магнитного поля Земли в данной работе использовались модели IGRF и EMM [Alken et al., 2021; The National Centers for Environmental Information, 2018]. Данные модели представляют магнитное поле Земли в виде разложения в ряд по сферическим гармоникам. Модель IGRF описывает главное магнитное поле Земли, создаваемое источниками в жидком ядре Земли, и использует 13 сферических гармоник для аппроксимации скалярного магнитного потенциала. Более подробная модель EMM описывает сумму главного и аномального магнитного поля, которое создается источниками в земной коре. Таким образом, вычитая из модуля поля модели EMM модуль поля модели IGRF, можно получить данные об аномальной составляющей полного магнитного поля

$\Delta T_1^A\left({\overrightarrow{r}}_i\right)=T^{EMM}\left({\overrightarrow{r}}_i\right)-T^{IGRF}\left({\overrightarrow{r}}_i\right)$, (3)

где $\overrightarrow{r}=(x,y,z)$ – координаты точки измерения, $T\left(\overrightarrow{r}\right)=\sqrt{H_X^{}{\left(\overrightarrow{r}\right)}^2+H_Y^{}{\left(\overrightarrow{r}\right)}^2+H_Z^{}{\left(\overrightarrow{r}\right)}^2}$ – модуль соответствующего магнитного поля, а число точек измерения i = 1…N_st. Поле (3) можно получить, например, если из измеренной протонным магнитометром сетки данных вычесть модуль главного магнитного поля модели IGRF.

С помощью простой итерационной процедуры по таким модульным данным можно рассчитать векторное аномальное магнитное поле, при условии, что компоненты главного поля, $H_X^{IGRF}\left({\overrightarrow{r}}_i\right),H_Y^{IGRF}\left({\overrightarrow{r}}_i\right),H_Z^{IGRF}\left({\overrightarrow{r}}_i\right)$, известны.

Пусть ${\overrightarrow{H}}_0$ и ${\overrightarrow{H}}_1$ есть главное и аномальное магнитное поле, соответственно. В нашем случае $T^{EMM}=\left|{\overrightarrow{H}}_0+{\overrightarrow{H}}_1\right|$, а $T^{IGRF}=\left|{\overrightarrow{H}}_0\right|$. Воспользуемся тем, что модуль главного магнитного поля, как правило, много больше модуля аномального поля [Яновский, 1978; Lourenco and Morrison, 1973]. Тогда имеем

$\begin{array}{c}\Delta T_1^A\left({\overrightarrow{r}}_i\right)=T^{EMM}-T^{IGRF}=\left|{\overrightarrow{H}}_0+{\overrightarrow{H}}_1\right|-\left|{\overrightarrow{H}}_0\right|\approx \\ \approx \left|{\overrightarrow{H}}_0\right|\sqrt{1+2{\overrightarrow{H}}_0{\overrightarrow{H}}_1/H_0^2}-\left|{\overrightarrow{H}}_0\right|\approx {\overrightarrow{H}}_0{\overrightarrow{H}}_1/\left|{\overrightarrow{H}}_0\right|.\end{array}$ (4)

Следовательно, в результате операции (4) получается компонента аномального поля в направлении глобального поля. Если глобальное поле вблизи полюса направлено практически вертикально, то в результате мы получаем с большой точностью вертикальную проекцию аномального магнитного поля, $H_Z^A\left({\overrightarrow{r}}_i\right)=\Delta T_1^A\left({\overrightarrow{r}}_i\right)$.

Далее, по данным $H_Z^A\left({\overrightarrow{r}}_i\right)$ с помощью эволюционного алгоритма рассчитывается компоненты поля в плоскости $H_X^A\left({\overrightarrow{r}}_i\right)$, $H_Y^A\left({\overrightarrow{r}}_i\right)$. Затем, зная все компоненты аномального поля,  $H_X^A\left({\overrightarrow{r}}_i\right)$, $H_Y^A\left({\overrightarrow{r}}_i\right)$, $H_Z^A\left({\overrightarrow{r}}_i\right)$, а также все компоненты главного поля, $H_X^{IGRF}\left({\overrightarrow{r}}_i\right)$, $H_Y^{IGRF}\left({\overrightarrow{r}}_i\right),H_Z^{IGRF}\left({\overrightarrow{r}}_i\right)$, находим модуль полного поля на первой итерации

${\overrightarrow{H}}^{Total}\left({\overrightarrow{r}}_i\right)={\overrightarrow{H}}^{IGRF}\left({\overrightarrow{r}}_i\right)+{\overrightarrow{H}}^A\left({\overrightarrow{r}}_i\right)$, (5)

$T^{Total}\left({\overrightarrow{r}}_i\right)=\sqrt{H_X^{Total}{\left({\overrightarrow{r}}_i\right)}^2+H_Y^{Total}{\left({\overrightarrow{r}}_i\right)}^2+H_Z^{Total}{\left({\overrightarrow{r}}_i\right)}^2}$. (6)

На второй итерации вычитаем из полного поля модели EMM модуль полного поля (6), найденного на первой итерации

$\Delta T_2^A\left({\overrightarrow{r}}_i\right)=T^{EMM}\left({\overrightarrow{r}}_i\right)-T^{Total}\left({\overrightarrow{r}}_i\right)$. (7)

Новая разница модулей полей (7) вновь подается на вход эволюционного алгоритма в качестве приближения к вертикальной компоненте искомого поля, после чего вышеописанная процедура (5)–(7) повторяется. Этот итерационный процесс продолжается до тех пор, пока разница между известным и расчетным полным магнитным полем не будет удовлетворять условию $\sum _i\left|T^{EMM}\left({\overrightarrow{r}}_i\right)-T^{Total}\left({\overrightarrow{r}}_i\right)\right|/N_{st}\le \delta$, где d − малая величина в сравнении с характерным значением аномального магнитного поля.

На первой стадии расчета проводилось определение оптимального количества диполей и их магнитных моментов по описанной в пункте 3 процедуре. Как показывает рис. 4, для восстановления векторного магнитного поля использовались N_d = 26 диполей с магнитным моментом M = 30∙10¹² А∙м². Размещение диполей проводилось в области с размерами L_x = 300 км, L_y = 300 км и L_z = 60 км, высота сканирования магнитного поля h = 100 м.

 

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

 

Источниками аномального магнитного поля в Земной коре обычно являются залежи, включающие в свой состав мелкие магнитные частицы, например, магнетит [Яновский, 1978]. Соответственно, намагниченность дипольных источников, используемых в алгоритме, по крайней мере, не должна превышать характерных значений намагниченности естественных горных пород.

В рассматриваемом случае намагниченность набора точечных диполей, распределенных в немагнитном объеме, можно оценить по формуле J = MN_d/V_a, где V_a = L_xL_yL_z − объем области размещения диполей под поверхностью земли. Так, набору диполей с параметрами M = 30∙10¹² А∙м² и N_d = 26, обозначенному окружностью на рис. 4, соответствует намагниченность J = 0.14 А/м, при этом характерное значение намагниченности горных пород составляет J_R = 3∙А/м [Montesinos et al., 2016]. Сравнивая J_R с полной намагниченностью дипольных источников J, находим отношение J/J_R = 0.046, что указывает на то, что намагниченность набора диполей J соответствует намагниченности горной породы с низкой концентрацией магнетита.

На рис. 5а показана область полуострова Ямал с координатами: 69° – 72° N и 66° – 74° E, для которой проводился расчет векторного поля по аномальной составляющей геомагнитного поля. Измерения проводились на сетке 23×23 точек, показанной на рис. 5б. Всего для расчета компонент поля проводилось усреднение по 100 запускам алгоритма. Итоговая величина невязки целевой функции f = 0.05.

 

Рис. 5. (а) − Исследуемая область; (б) − исходные модульные данные аномального магнитного поля DT, полученные по уравнению (4) в исследованной области. (в, г, д) – найденные с помощью вычислений компоненты H_x, H_y, H_z магнитного поля; (е, ж, з) − истинные H_x, H_y, H_z компоненты магнитного поля, полученные из разности моделей EMM и IGRF, ${\mathit{H}}_{\mathbf{i}}\mathbf{=}{\mathit{H}}_{\mathbf{i}}^{\mathbf{E}\mathbf{M}\mathbf{M}}\mathbf{-}{\mathit{H}}_{\mathbf{i}}^{\mathbf{I}\mathbf{G}\mathbf{R}\mathbf{F}}$, i = x, y, z.

 

Рисунок 5б показывает, что аномальная составляющая геомагнитного поля DT, вычисленная по уравнению (4) в исследуемой области, принимает значения от – 230 нТ до 284 нТ. С помощью эволюционного алгоритма по данным на рис. 5б было рассчитано векторное магнитное поле, компоненты которого показывают рис. 5в–д. Видно, что контуры распределения вертикальной компоненты H_z, показанные на рис. 5д, совпадают контурами DT на рис. 5б, что находится в согласии с приближением (4). Далее, на рис. 5в–д показаны истинные распределения компонент аномального магнитного поля, полученные напрямую из разности известных компонент моделей EMM и IGRF, $H_i=H_i^{EMM}-H_i^{IGRF}$, i = x, y, z. Как показывают рис. 5в–д, контуры компонент рассчитанного аномального магнитного поля хорошо совпадают с истинными компонентами поля, рис. 5е–з.

Интересно проследить показанную на рис. 6 эволюцию разности $T^{EMM}\left({\overrightarrow{r}}_i\right)-T^{Total}\left({\overrightarrow{r}}_i\right)$ между модулем поля EMM и модулем поля, восстанавливаемого с помощью эволюционного алгоритма, в зависимости от номера итерации.

 

Рис. 6. Эволюция разности (в нТ) между истинным и рассчитанным с помощью эволюционного алгоритма распределением магнитного поля, ${\mathit{T}}^{\mathbf{E}\mathbf{M}\mathbf{M}}\mathbf{\left(}{\overrightarrow{\mathbf{r}}}_{\mathbf{i}}\mathbf{\right)}\mathbf{-}{\mathit{T}}^{\mathbf{T}\mathbf{o}\mathbf{t}\mathbf{a}\mathbf{l}}\mathbf{\left(}{\overrightarrow{\mathbf{r}}}_{\mathbf{i}}\mathbf{\right)}$, (а) − на первой итерации, (б) − на второй итерации и (в) − на третьей итерации по уравнениям (5) – (7).

 

Как показывает рис. 6а, на первой итерации максимальная разность достигает значения в 140 нТ, на второй итерации, рис. 6б, разность уменьшается до значений порядка 70 нТ, но уже на третьей итерации, рис. 6в, максимальная разность $T^{EMM}\left({\overrightarrow{r}}_i\right)-T^{Total}\left({\overrightarrow{r}}_i\right)$ не превышает 30 нТ. Можно показать, что с увеличением числа итераций в процедуре (5) – (7) невязка $\sum _i\left|T^{EMM}\left({\overrightarrow{r}}_i\right)-T^{Total}\left({\overrightarrow{r}}_i\right)\right|/N_{st}$ стремится к нулю, что позволяет увеличить точность восстановления векторного магнитного поля.

5. ЗАКЛЮЧЕНИЕ

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

Тестирование алгоритма на синтетических данных показало, что алгоритм эффективен в построении карт распределений полного векторного магнитного поля на некоторой поверхности.

В работе показано, что с помощью предлагаемой итерационной процедуры можно рассчитать векторное аномальное магнитное поле по измеренным данным для разности модулей полей DT в некоторой области. Тестовые расчеты по данным моделей IGRF и EMM показали эффективность разработанного метода, а определенная итерационной процедурой разница между известным и расчетным магнитными полями становится ниже 30 нТл уже на третьей итерации эволюционного алгоритма. На данный момент ведутся работы по тестированию алгоритма на базе измеренных локальных аномальных магнитных полей и дальнейшему применению разработанного метода для создания локальных моделей геомагнитного поля.

ФИНАНСИРОВАНИЕ

Исследования выполнены в рамках государственного задания института.

Об авторах

Р. А. Рытов

Институт земного магнетизма, ионосферы и распространения радиоволн им. Н.В.Пушкова РАН (ИЗМИРАН)

Автор, ответственный за переписку.
Email: ruslan.rytov2017@ya.ru
Россия, Москва, Троицк

Н. А. Усов

Институт земного магнетизма, ионосферы и распространения радиоволн им. Н.В.Пушкова РАН (ИЗМИРАН)

Email: usov@obninsk.ru
Россия, Москва, Троицк

В. Г. Петров

Институт земного магнетизма, ионосферы и распространения радиоволн им. Н.В.Пушкова РАН (ИЗМИРАН)

Email: vgpetrov2018@mail.ru
Россия, Москва, Троицк

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

Дополнительные файлы