Физика Земли, 2023, № 2, стр. 20-35

ВВЕДЕНИЕ

За последние 10–20 лет объем информации о физических полях и топографии Земли увеличился на несколько порядков. Благодаря различным спутниковым миссиям непрерывным потоком поступают данные о планетах Солнечной системы, при этом орбиты космических аппаратов фактически “заметают” некоторые области трехмерного пространства, обеспечивая, таким образом, достаточно высокую степень разрешения съемки. Как уже было отмечено в нашей предыдущей работе [Степанова и др., 2020], материалы дистанционного зондирования Земли и других небесных тел используются при уточнении теории внутреннего строения нашей планеты и ей подобных, незаменимы при прогнозировании глобального изменения климата и изучении разнообразных геофизических процессов.

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

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

Такой тезис, на первый взгляд, кажется парадоксальным, поскольку локальные аномалии гравитационного (магнитного, теплового и т.п.) поля, как правило, “индивидуальны”: у каждой имеется свой собственный, отличный от других, “портрет”. Но можно объединить отдельные аномалии в группы и подгруппы и назвать такие множества “классами” и “подклассами” объектов, а операции над такими объектами отнести к соответствующим классам отображений – подобно тому, как это делается в математической теории категорий и в программных средах (языках) объектного программирования. К такому выводу авторы настоящей работы пришли на основании многолетнего опыта в решении различных линейных и нелинейных обратных и некорректных задач математической геофизики.

При интерпретации геофизических данных применяются различные методы и подходы (см., например, [Portniaguine, 1999; 2002; Arkani-Hamed, 2002; Acuna, 1999; Михайлов и др., 2019; Wang et al., 2019; Долгаль, 2015; 2017; Uieda et al., 2015; Schattner et al., 2019]).

Восстановление источников сигнала по информации об элементе поля и его градиенте [Wang et al., 2020a; 2020b] представляется нам особенно перспективным.

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

В работах авторов, посвященных применению метода линейных интегральных представлений для решения некорректных обратных задач геофизики, ранее не рассматривались постановки для уравнений второго порядка параболического и гиперболического типов. Между тем, возможности указанного метода не ограничиваются уравнениями Лапласа и Пуассона, которым, как известно, удовлетворяют компоненты гравитационного, электрического, магнитного и т.д. полей при определенных условиях. Для того, чтобы восполнить этот “пробел”, мы попытались в настоящей работе предложить аппроксимационную конструкцию для математического моделирования магнитного поля Земли, в которой в явном виде будет учитываться время. Необходимо подчеркнуть, что вне источников поля (т.е. вне внешнего жидкого ядра Земли и других космических объектов, обладающих внутренним генератором магнитного поля) методика разложения сигналов в ряд по сферическим гармоникам, представления в виде интегралов Фурье, вейвлетов и т.п. [Фрик и др., 2020; Титов и др., 2020; Казанцев, 2019] вполне допустимы и могут хорошо “работать” при решении задач анализа и синтеза. Но адекватное реальности моделирование магнитного поля Земли невозможно получить без учета магнитогидродинамических соотношений [Reshetnyak, 2015]. Поэтому актуальной, на наш взгляд, является проблема нахождения устойчивых приближенных решений системы уравнений магнитной гидродинамики, в которой для магнитного поля могут считаться справедливыми некоторые интегральные представления. Мы, таким образом, получаем некоторую систему обыкновенных дифференциальных уравнений, которой удовлетворяет поле скоростей космических частиц – солнечного ветра. Аналогичным образом можно описывать потоки заряженной жидкости (плазмы) и внутри планет. Первым шагом на пути такого моделирования полей является учет времени в качестве независимой переменной (или параметра модели, поскольку при построении решений уравнений в частных производных продуктивным является метод спуска по одной из переменных, от которых эти решения зависят [Владимиров, 1981]).

В наших предыдущих работах мы показали, что глобальный вариант метода линейных интегральных представлений в рамках структурно-параметрического подхода может применяться для построения аналитического продолжения магнитного поля Марса вниз, в сторону источников [Gudkova et al., 2020; 2021; Salnikov et al., 2021]. На Марсе динамо в настоящее время отсутствует, и спутниковые данные свидетельствуют о фрагментарной картине остаточной намагниченности на Красной планете. Как мы уже отмечали в приведенных выше работах, причина исчезновения магнитного поля Марса не установлена.

В настоящей работе мы предприняли попытку аппроксимации аномального магнитного поля Земли в районе Японского моря, согласно данным модели EMAG2v3 [Meyer et al., 2017]. Модель EMAG3v3 получена объединением низкочастотных спутниковых данных, в основном представленных сферической моделью MF7, вычисленной по измерениям спутниковой миссии CHAMP в период с 2007 по 2010 гг., с большим объемом инструментальных морских и аэромагнитных съемок. Пространственное разрешение модели составляет 2', и мы выбирали именно такое значение шага дискретизации при постановке математического эксперимента.

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

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

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

ОПИСАНИЕ МЕТОДИКИ

Мы неоднократно подчеркивали в предыдущих работах [Stepanova, 2008; Раевский и др., 2015а; 2015б], что локальные интегральные представления элементов физических полей и других сигналов строятся для небольших по площади территорий: не более 10⁴ км², в отличие от глобальных (линейные аналитические аппроксимации строятся для всей Земли) или региональных (линейные аналитические аппроксимации строятся для больших по площади территорий – порядка (10⁵–10⁶) км²) [Stepanova, 2009].

В данной работе будут рассматриваться зависящие от времени физические поля и другие сигналы в декартовой системе координат – можно полагать, что мы ограничиваемся пока только локальным вариантом ввиду сложности задачи аналитического описания нестационарного поля в глобальном масштабе. Однако такой подход не ограничивает общности: в принципе, можно решать и глобальные задачи в общеземной системе координат. Но тогда эквивалентные по внешнему полю источники окажутся принадлежащими слишком абстрактному классу объектов. С математической точки зрения это – упрощение, некая идеализация, нулевое приближение к решению нелинейной обратной задачи [Ягола и др., 2014]. И такое решение имеет право на существование, поскольку при изучении планет строятся, как правило, метрологические аппроксимации полей: количественная информация об источниках совершенно не доступна [Страхов, Степанова, 2002]. В работе [Степанова и др., 2019] мы предложили применять для интерпретации данных гравимагниторазведки метод линейных интегральных представлений повышенной размерности. При таком подходе, как показал математический эксперимент, особенности векторных полей выявляются с большей степенью точности при построении метрологических аппроксимаций.

Если эллиптичность не учитывается, то планета считается шаром радиуса ${{R}_{0}}$.

В том случае, когда эллиптичность нужно учитывать, Земля представляется в виде эллипсоида с полуосями a, b, c.

Если известны геодезические координаты точки наблюдения, то декартовы можно выразить через геодезические по известным формулам [Шимбирев, 1975]:

где: X, Y, Z – пространственные прямоугольные координаты точки, м; B, H, L – геодезические координаты (широта и долгота, рад., высота, м); N – радиус кривизны первого вертикала, м; e – эксцентриситет эллипсоида.

Радиус кривизны первого вертикала и квадрата эксцентриситета эллипсоида вычисляются по формулам:

где: a − большая полуось эллипсоида, м; α − сжатие эллипсоида.

“Идеализированную” Землю будем считать внутренностью эллипсоида с полуосями a, b, c, а реальную Землю представим как область, ограниченную кусочно-непрерывной замкнутой поверхностью S, мало отклоняющейся от этого эллипсоида и содержащей его внутри себя. Предположим, что в некоторой (произвольной) совокупности точек ${{x}^{{\left( i \right)}}},i = 1,2,...,N,$ на поверхности S заданы приближенные значения гармонической (во внешности эллипсоида) функции $G(x)$:

Здесь $G\left( x \right)$ имеет следующее интегральное представление:

Относительно функций $\sigma $ и $w$ ставится следующая условно вариационная задача:

Функции $p_{1}^{2}(\vartheta ,\varphi ),\;p_{2}^{2}(\vartheta ,\varphi )$ в вариационной постановке (5) должны отражать свойства области задания эквивалентных по внешнему полю носителей масс. В частности, эти функции могут принимать заданные значения на границе области, внутри или вне которой нужно найти распределение элементов поля. В то же время, сами указанные функции могут быть неизвестными, и их определение можно сравнить с механизмом адаптивной регуляризации: элементы матрицы системы линейных алгебраических уравнений (СЛАУ), к решению которой редуцируется вариационная постановка (5), зависят от $p_{1}^{2}(\vartheta ,\varphi ),\;p_{2}^{2}(\vartheta ,\varphi )$ . Можно рассмотреть, таким образом, некоторое семейство вариационных задач по определению как неизвестных распределений масс, так и областей задания источников. Подобные вариационные постановки относятся к нелинейным обратным задачам, что существенно осложняет процесс поиска решения. Наша методика построения аналитических моделей зависящего от времени магнитного поля Земли может быть описана следующим образом.

Решение вариационной задачи (5) при известных функциях $p_{1}^{2}(\vartheta ,\varphi ),\;p_{2}^{2}(\vartheta ,\varphi )$ эквивалентно решению системы линейных алгебраических уравнений [Страхов, Степанова, 2002] :

где вектор ${{f}_{\delta }}$ имеет компоненты ${{f}_{{i,\delta }}}$, а матрица ${\text{A}} = {{{\text{A}}}^{{\text{T}}}} \geqslant 0$ имеет элементы

Для $\sigma $ и $w$ получаем представление:

где положено

Если функции $p_{1}^{2}(\vartheta ,\varphi ),\;p_{2}^{2}(\vartheta ,\varphi )$ тождественно на сфере равны 1, то для элементов матрицы СЛАУ получаем следующие выражения:

где $Q_{i}^{{(1)}}(\xi )$ и $Q_{i}^{{(2)}}(\xi )$ определяются выражениями:

Здесь ${{h}_{i}} = \frac{{{{R}_{0}}}}{{{{r}_{i}}}},\,\,{{h}_{j}} = \frac{{{{R}_{0}}}}{{{{r}_{j}}}}$ , а ${{\alpha }_{{ij}}}$ – это угол между векторами ${{x}_{i}}$ и ${{x}_{j}}$. Выражения (10) были получены в предположении, что вектор ${{x}_{i}}$ параллелен оси Oz, в то время как вектор ${{x}_{j}}$ лежит в плоскости zOx (в случае сферы мы всегда можем выбрать такую систему координат, в отличие от эллипсоида). Тогда радиус-векторы i-ой и j-ой точек будут иметь координаты ${{x}_{i}} = ({{r}_{i}},0,0)$ и ${{x}_{j}} = ({{r}_{j}}\sin {{\alpha }_{{ij}}},\,\,0,{{r}_{j}}\cos {{\alpha }_{{ij}}})$.

Функция $F(\varphi ,k)$ в (10) представляет собой эллиптический интеграл первого рода.

Если вместо сферы поверхностью интегрирования является эллипсоид с полуосями a, b, c (т.е. верно интегральное представление (4) и рассматривается соответствующая вариационная постановка), то следует перейти в пространство других декартовых координат.

А именно, положим:

Эллипсоид с указанными полуосями переходит при такой замене координат в сферу единичного радиуса, а точки наблюдения – в точки штрихованного пространства по формулам (12), внешние по отношению к единичной сфере.

Преобразование (12) не является ортогональным и, по этой причине, расстояния между точками не сохраняются. Если решить вариационную задачу для “штрихованной” единичной сферы, а затем перейти в “нормальное” декартово пространство, то элементы матрицы (10) принимают вид:

Здесь ${{\tilde {h}}_{i}} = \frac{1}{{r_{i}^{'}}},$ $\tilde {h}_{j}^{'} = \frac{1}{{r_{j}^{'}}},$ $r_{i}^{'} = \sqrt {x_{{1,i}}^{{'2}} + x_{{2,i}}^{{'2}} + x_{{3,i}}^{{'2}}} ,$ $r_{j}^{'} = \sqrt {x_{{1,j}}^{{'2}} + x_{{2,j}}^{{'2}} + x_{{3,j}}^{{'2}}} $, а ${{\tilde {\alpha }}_{{ij}}}$ – это угол между векторами $x_{i}^{'}$ и $x_{j}^{'}$.

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

Безусловно, при таком подходе невозможно учесть все топологические особенности планеты, но представление элементов матрицы СЛАУ в аналитическом виде существенно упрощает процесс нахождения устойчивого приближенного решения обратной задачи. Как видно из выражений (12) и (13), время в таких моделях никак не фигурирует.

При задании пространственных координат точек наблюдения необходимо руководствоваться формулами (1)–(2), связывающими геодезические (или сферические) и декартовы координаты точки наблюдения в общеземной системе координат.

Нестационарность магнитного поля может учитываться по-разному: все определяется конкретной постановкой. В последние 10–20 лет для интерпретации данных о магнитных полях Земли и планет особенно популярным среди зарубежных исследователей стал так называемый метод сингулярных магнитных диполей (ESD) [Whaler et al., 2005; Langlais et al., 2004]. Существует также вариант указанного метода с зависящими от времени компонентами вектора намагниченности [Oliveira et al., 2015]. При этом авторы выбирали наиболее простую форму зависимости – линейную.

В общем случае, поведение магнитного поля планеты (как внешнего, так и внутреннего) определяется из решения достаточно сложной системы уравнений магнитной гидродинамики [Арнольд, Хесин, 2007]:

что является так называемым линеаризованным и укороченным уравнением Эйлера.

В (14): $\bar {v}(\bar {r},t),\;\bar {B}(\bar {r},t)$ - соответственно поле скоростей и магнитная индукция как функции радиус-вектора и времени; $\nabla p$ – градиент давления в среде; $\Delta $ – оператор Лапласа. Квадратные скобки означают векторное произведение двух векторов, а фигурные скобки – это скобка Пуассона двух векторных полей, определяемая как:

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

Если мы интерпретируем данные о магнитном поле планеты, полученные в результате дистанционного зондирования, то система (14) (или подобная ей) должна учитывать взаимодействие внутреннего (достаточно сильного) поля с магнитным полем, создаваемым потоком заряженных частиц (солнечного ветра), если же нас интересует сам процесс генерации динамо, то поле скоростей в (14) – это поток раскаленной жидкости в ядре.

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

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

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

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

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

Две составляющие аналитической модели магнитного поля могут показаться избыточными, на первый взгляд. Однако магнитное поле, как мы уже подчеркивали, не может, в общем случае, рассматриваться как решение уравнения Лапласа (в области, не занятой магнитными массами) или Пуассона (в коре планеты, но не в ядре при наличии магнитного динамо).

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

Опишем теперь, как строится вторая “компонента” математической модели магнитного поля Земли. В настоящей работе мы полагаем, что время – это некий непрерывно меняющийся от нуля до бесконечности параметр модели . Параметр ${{a}^{2}}$ можно считать достаточно большим: если ${{a}^{2}} \to \infty $, то решение приведенного ниже уравнения параболического типа будет стремиться к решению уравнения Лапласа. Ну, а в случае системы (14) указанный параметр связан с коэффициентом магнитной вязкости (магнитным числом Рейнольдса), но мы пока не рассматриваем решение системы (14) в целом.

Как известно, метод спуска по одной из независимых переменных [Владимиров, 1981] позволяет получить решение краевой задачи в пространстве N переменных, если найдено соответствующее решение в N + 1-ом пространстве. Поэтому переход от вектор-функции, зависящей от трех пространственных переменных, к вектор-функции в четырехмерном пространстве (четвертая координата – это время) нам кажется оправданным при построении метрологических аппроксимаций элементов физических полей и топографии.

Если считать Землю (или другую планету) нижним полупространством, ограниченным плоскостью ${{x}_{3}} = 0$ в декартовой системе координат (${{x}_{1}}~{{x}_{{2,}}}~{{x}_{3}}$), то решение первой краевой задачи для уравнения теплопроводности (или диффузии, важно – уравнения параболического типа с постоянными коэффициентами в верхнем полупространстве) можно записать в следующем виде:

Формула (15) получена с помощью нечетного продолжения фундаментального решения уравнения теплопроводности во всем пространстве на нижнее полупространство (см., например, [Будак и др., 1980]). Решение неоднородного уравнения теплопроводности во всем пространстве имеет вид:

Представления решения краевой задачи (15) или уравнения относительно вектор-функции $\bar {B}({\mathbf{x}},t)$ (16) позволяют нам применить метод линейных интегральных представлений [Страхов, Степанова, 2002] в локальном варианте, когда Земля представляется в виде нижнего полупространства, а компоненты магнитной индукции известны в верхнем полупространстве.

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

Обратимся к формуле (15) и попробуем сформулировать задачу по поиску минимума некоторого функционала, опираясь на соотношения (5). Вариационные постановки, возникающие при применении различных версий метода линейных интегральных представлений, предполагают наличие у минимизируемой по норме ${{L}_{2}}({{M}_{r}})$ функции (функций, вектор-функций и т.п.) определенных свойств гладкости. Здесь ${{L}_{2}}({{M}_{r}})$ – гильбертово пространство интегрируемых с квадратом функций в областях задания источников поля ${{M}_{r}},\;r = 1,...,R.$

Если вектора магнитной индукции удовлетворяют уравнениям параболического типа, то интеграл по времени от произведения мощности источника и функции Грина для всего пространства берется от нуля до того момента времени, в который было осуществлено измерение сигнала (t, например). Для уравнения Лапласа мы всегда ставили вариационные задачи при фиксированном наборе множеств ${{M}_{r}},~r = 1, \ldots ,R.$ Мы предлагаем определять неизвестную мощность источников из решения следующей вариационной постановки:

Решением вариационной задачи (17) будет служить функция:

Для элементов матрицы, соответствующей СЛАУ, получим выражения:

где $\Phi (\zeta ) - $ функция ошибок:

При решении вариационных задач, основанных на соотношениях (15), можно считать, что функция $\frac{{\partial {{V}_{a}}\left( x \right)}}{{\partial {{x}_{3}}}}$, гармоническая в полупространстве x₃ > –H (или ${{x}_{3}} < H$ ) и связанная с потенциалом аномального магнитного поля ${{V}_{a}}(x)$, имеет следующее интегральное представление:

Здесь интегрирование выполняется по плоскости, на которой в данный момент времени находятся источники нестационарного поля (${{x}_{3}} = - H$).

Представим неизвестную в (20) функцию в следующем виде:

В этом случае задача о нахождении неизвестной функции $\psi (\tau ,{{x}_{1}},{{x}_{2}})$ сводится к определению вектора $\lambda = ({{\lambda }_{1}}, \ldots {{\lambda }_{N}})$.

Мы получаем характерную для постановок в рамках метода линейных интегральных представлений СЛАУ (систему линейных алгебраических уравнений):

в которой А есть (N × N) – матрица со свойством и элементами ${{a}_{{pq}}},\,\,1 \leqslant p,\,\,q \geqslant N$: а ${{f}_{\delta }}$ есть N – вектор с компонентами ${{f}_{{i,\delta }}}$ , определенными по (3.9); λ – есть N – вектор с компонентами ${{\lambda }_{i}}$.

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

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

1) По заданному набору точек наблюдения в трехмерном пространстве и заданным начальным приближениям к неизвестным весовым функциям ${{p}_{r}},{\kern 1pt} r = 1,2$, строится модифицированная S-аппроксимация в региональном варианте. Декартовы координаты при этом связаны с геодезическими соотношениями (1)–(2).

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

3) Весовые функции меняются в соответствии с определенным исследователем правилом, процесс повторяется до тех пор, пока не будет достигнута необходимая точность.

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

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

При проведении тестовых расчетов был выбран следующий эмпирический закон изменения весовых функций:

Здесь: $r = \sqrt {{{{\left( {\xi - {{x}_{0}}} \right)}}^{2}} + {{{\left( {\eta - {{y}_{0}}} \right)}}^{2}} + {{{\left( {\varsigma - {{z}_{0}}} \right)}}^{2}}} ,$ $P = \left( {\xi ,\eta ,\varsigma } \right)$ – текущая точка на поверхности, по которой производится интегрирование; $M = \left( {{{x}_{0}},{{y}_{0}},{{z}_{0}}} \right)$ – координата центра масс рассматриваемого набора данных. Показатели степени m и n предполагаются равными целым положительным числам.

При нахождении S-аппроксимаций, иначе говоря – функций $\sigma \left( \xi \right)$ и $w\left( \xi \right)$, основная вычислительная проблема состоит в нахождении устойчивого приближенного решения (согласованного с имеющейся априорной информацией о векторах f и $\delta f$) системы линейных алгебраических уравнений (7)–(8).

Для решения систем (6), (7), (13) и (22)–(24) мы применили два метода, предложенные авторами в работах [Stepanova et al., 2020; Pan et al., 2020; Gudkova et al., 2020; Salnikov et al., 2020; 2021]. Первый из методов решения, плохо обусловленных СЛАУ, мы будем обозначать как БМХР (блочный метод регуляризации Холецкого), а второй – как УБМ (усовершенствованный блочный метод). В основе указанных методов лежит механизм адаптивной регуляризации. Сначала решается система уравнений, матрица которой строится по всему набору данных, а затем некоторые блоки матрицы “выключаются” из вычислительного процесса, если соответствующие компоненты вектора решения не оказывают значимого влияния на результат в целом. Мы протестировали оба метода на большом числе модельных и практических примеров, созданных по материалам гравимагниторазведки. В основном, при интерпретации гравиметрических и магнитометрических данных приходится решать СЛАУ, имеющих матрицы с диагональным преобладанием. Блоки матрицы, в которых содержится информация о взаимодействии далеко друг от друга расположенных носителей масс, отличаются существенно меньшим, по сравнению с диагональными, рангом. Поэтому можно сделать однозначный вывод о нецелесообразности выбора одного и того же параметра регуляризации для диагональных и внедиагональных блоков плохо обусловленной системы линейных алгебраических уравнений. В БМХР для каждого из блоков матрицы СЛАУ выполняется регуляризация разложения Холецкого, а в УБМ применяются как итерационные алгоритмы для поиска решения, так и прямые (основанные на регуляризации LU-разложения несимметричной матрицы). Усовершенствованный блочный метод позволяет решать систему и с единственным блоком. Автоматически выбор количества блоков матрицы и вектора параметров регуляризации пока не осуществляется ввиду высокой степени неустойчивости приближенных решений обратных линейных (а тем более, нелинейных) обратных задач геофизики. На данный момент возможны исключительно метрологические аппроксимации гравитационного и магнитного полей, поскольку количественная информация об источниках полей отсутствует.

АПРОБАЦИЯ АЛГОРИТМОВ И ПРОГРАММ НОВОЙ КОМБИНИРОВАННОЙ МЕТОДИКИ ПОСТРОЕНИЯ АНАЛИТИЧЕСКИХ МОДЕЛЕЙ МАГНИТНОГО ПОЛЯ ПО ДАННЫМ МОДЕЛИ EMAG3V3 В РАЙОНЕ ЯПОНСКОГО МОРЯ

В работе [Родников и др., 2014] приводится сравнительно подробное описание морфологических особенностей магнитного поля в районе Японского моря. В юго-западной части Японского моря аномалии имеют субширотный характер и северо-западное и субмеридиональное простирание. В южной части моря наблюдается слабоинтенсивное знакопеременное поле, при этом аномалии магнитного поля достигают 300 нТл. В юго-западной части интенсивность аномалий не превосходит в среднем 100 нТл, иногда встречаются значения 300 нТл. Восточная часть Японского моря отличается слабоинтенсивными аномалиями до ±50 нТл, однако отрицательные локальные аномалии порой достигают значений 200 нТл. Как указывают авторы [Родников и др., 2014], существенного различия в структуре поля районов, непосредственно примыкающих к северной части о. Хонсю и удаленных от него, не отмечается. Здесь встречаются высокочастотные полосовые аномалии субширотного и северо-восточного простираний, амплитуда которых составляет ≥400 нТл. Положительным аномалиям можно поставить в соответствие массивы неогеновых и четвертичных вулканических пород. На северной оконечности острова преобладают отрицательные аномалии интенсивностью –400 нТл преимущественно меридионального простирания. Для примыкающих к северной части острова Хонсю участках северо-западной части Тихого океана можно отметить отрицательные (–400 нТл) и положительные (+200 нТл) линейные аномалии северо-восточного простирания.

Карта аномалий магнитного поля в районе Японского моря изображена на рис. 1 (согласно модели EMAG3v3).

Мы выбрали для тестирования новой методики полигон, размеры которого составили 12° по широте и 12° по долготе. Полигон был разбит на три подобласти по широте с целью выявления “сильных” и “слабых” мест нового алгоритма интерпретации магнитных данных. Для каждой из трех подобластей были сформированы по два массива точек наблюдения: один массив включал в себя значения географических координат и значения модуля аномального магнитного поля на уровне моря согласно модели EMAG3v3, а второй – отсортированные по значениям модуля магнитной индукции точки). В таблице 1 первые три строки соответствуют массивам первого типа, а последние три строки – второго типа. На первом этапе построения математической модели магнитного поля мы применили региональный вариант модифицированного метода S-аппроксимаций с указанными в предыдущем разделе статьи весовыми функциями.

Магнитное поле аппроксимировалось суммой простого и двойного слоев, распределенных на поверхности эллипсоида, охватываемого земным шаром, (радиус Земли полагался равным 6378 км, полуоси эллипсоида – 6358–6370 км).

Результаты решения СЛАУ приведены в табл. 1; в последнем столбце этой таблицы содержатся значения показателя качества решения (фактически – относительной погрешности) и время, за которое была достигнута заранее заданная точность решения.

Весовые функции рассматривались трех типов:

Наилучший результат был достигнут для второго типа, и в табл. 1 приводятся показатели качества решения именно для этого случая.

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

В таблице 1 показаны результаты первого этапа алгоритма для всех наборов тестовых данных. На рисунках мы решили показать аппроксимированное и исходное поле для некоторых примеров, поскольку прослеживается одна и та же тенденция: чем больше разброс координат в тестовом наборе, тем выше качество аппроксимации даже при четырехкратном увеличении числа неизвестных (с 5000 неизвестных до 20 000).

Шаг дискретизации в модели EMAG3v3 равен 2 мин, при выполнении расчетов мы задавали такой же.

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

На рис. 2а приводится карта изолиний магнитного поля в северной подобласти исследуемого региона, N = 10 000. На рис. 2б приводится карта изолиний поля, построенного в результате S-аппроксимации при выборе весовых функций (25) II-го типа.

На рис. 3а показано поле южной подобласти. Число точек в наборе полагалось равным 10 000.

Карта разности между полученными в результате аналитического продолжения на высоту 4 км и теоретическими значениями поля приведена на рис. 3б. Весовые функции III типа.

На рис. 4 приводится карта изолиний аппроксимированного магнитного поля с учетом эллиптичности при N = 20 000 (здесь N – число точек наблюдений, по которым строится аппроксимация) для южной подобласти, массив данных был отсортирован по возрастанию значений модуля аномального магнитного поля. Весовые функции II типа.

После того, как для каждой из подобластей была построена модифицированная S-аппроксимация магнитного поля на уровне моря и на высоте 4 км, мы перешли ко второму этапу алгоритма. Мы поставили задачу (21)–(24) для каждого из описанных выше наборов данных о поле. Время можно ввести в функцию, зависящую лишь от пространственных координат, следующим образом (безусловно, предлагаемая далее конструкция не должна ограничивать общность: если известны компоненты вектора магнитной индукции в реальные моменты времени, то время будет независимой переменной в системе дифференциальных уравнений, заданной в некоторой области пространства и на некотором временном интервале). Будем считать, что воображаемое время связано с координатами точки наблюдения следующим образом: ${{t}_{i}} = \frac{{{{{\left| {{{x}^{{\left( i \right)}}}} \right|}}^{2}}}}{{{{a}^{2}}}}$. В числителе здесь подразумевается квадрат расстояния от точки наблюдения до начала координат. Мы полагали ${{a}^{2}}$ равным 1000, однако такой подход не ограничивает общности. Если бы вместо уравнения параболического типа фигурировало волновое уравнение, то ${{a}^{2}}$ имело бы физический смысл квадрата скорости. В нашем представлении – это величина, характеризующая гипотетическую теплопроводность (скорее, рассеяние магнитного поля благодаря наличию вязкости). Подчеркнем, что на данном этапе разработки новой математической модели магнитного поля $~{{a}^{2}}$ рассматривается как параметр; обратная к $~{{a}^{2}}$ величина может считаться малым параметром, нулевое значение этого малого параметра соответствует “классической” постановке для уравнения Лапласа в области, не занятой магнитными массами. В принципе, каждой точке на регулярной или нерегулярной сети наблюдений можно сопоставить временной параметр произвольным образом. Если у всех точек временные координаты одинаковы, это может означать, что задано некоторое начальное распределение векторного поля. А впоследствии происходит процесс диффузии (или переноса гипотетического тепла) в пространстве и во времени. При генерации магнитного поля в недрах планеты или при взаимодействии внутреннего поля планеты с потоком заряженных частиц зависимость полей от времени определяется уравнением Эйлера, системой уравнений Максвелла и уравнением состояния среды. Мы предлагаем такую модель поля именно исходя из вида уравнений магнитной гидродинамики. Слишком простая, линейная, зависимость вектора намагниченности от времени могла рассматриваться ранее, но в настоящее время возникли новые вызовы, обусловленные стремительным ростом объемов получаемой в результате дистанционного зондирования Земли и других планет Солнечной системы информации.

Значения магнитного поля считались известными в сети точек на уровне моря и на высоте 4 км. По этим значениям мы определяли функцию ${{\Psi }}\left( {{{\tau }},{{{{\xi }}}_{1}},{{{{\xi }}}_{2}}} \right).$

На рис. 5 показано поле, аппроксимированное на втором этапе поля южной подобласти; N = 10 000; массив данных отсортирован. Параметр, играющий роль времени, был искусственно введен в постановку. Значения “времени” в i-ой точке определялись как: ${{t}_{i}} = \frac{{{{{\left| {{{x}^{{\left( i \right)}}}} \right|}}^{2}}}}{{{{a}^{2}}}}$. Параметр ${{a}^{2}}~$ полагался равным 1000.

В табл. 2 приводятся результаты аппроксимации магнитного поля для трех подобластей в районе Японского моря с помощью зависящего от времени варианта локальных модифицированных S-аппроксимаций.

Результаты расчетов на тестовых примерах показывают, что оба этапа комбинированной методики позволяют сохранить “естественный” характер магнитного поля: не наблюдается слишком резкого роста или убывания значений модуля аномального магнитного поля при изменении высоты точек наблюдения, т.е. поле не “распадается”.

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

ЗАКЛЮЧЕНИЕ

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