On Constructing Magnetic and Gravity Images of Mercury from Satellite Data

Толық мәтін

ВВЕДЕНИЕ

Информация об объектах Солнечной системы, полученная благодаря межпланетным миссиям, таким как Mariner-10, Messenger, BepiColombo, позволяет уточнять внутреннее строение планет и доказывать правоту выдвинутых ранее гипотез о происхождении физических полей, например, о внутреннем происхождении магнитного поля Меркурия. Благодаря данным Mariner-10 исследователи сделали вывод о внутреннем происхождении магнитного поля Меркурия [Ness et al., 1974; 1975; Johnson et al., 2012; Alexeev et al., 2010; Wicht et al., 2014 ]. Меркурий, как выяснилось, обладает целым рядом особенностей – в частности, достаточно сильным магнитным полем.

Миссия Messenger позволила изучить поверхность Меркурия, историю геологического развития, химический состав и магнитосферу [Hauck et al., 2007; Benkhoff et al., 2010; Anderson et al., 2007; 2008, Milillo et al., 2020]. Роль этого космического зонда в исследованиях Меркурия трудно переоценить: данные, собранные Messenger, свидетельствуют о существовании жидкой части ядра планеты [Anderson et al., 2010; 2011; 2012; 2013], что противоречит существовавшей до 1970-х годов гипотезе о наличии в ядре Меркурия лишь твердой фазы [Plageman, 1965]. Токи во внешнем, жидком, ядре создают внутреннее магнитное поле Меркурия. В работе [Anderson et al., 2012] было установлено, что магнитный экватор расположен в северной полусфере планеты, и этот факт послужил основанием для моделирования внутреннего поля Меркурия как поля единичного диполя. Значение дипольного момента Меркурия приведено в работе [Anderson et al., 2012]. Оказалось, что угол наклона диполя по отношению к оси вращения планеты невелик: он составляет около 0.8°. Изучались также и временные вариации магнитного поля Меркурия [Philpott et al., 2014].

Магнитное поле ближайшей к Солнцу планеты отличается весьма сложной структурой: есть составляющая, порождаемая движением горячей проводящей жидкости внутри ядра, а также поле плазмы в окружающем Меркурий пространстве. Однако по отдельности эти “компоненты” рассматривать нельзя: поле ядра активно взаимодействует с солнечным ветром. В магнитосфере Меркурия имеются области пересоединения магнитного поля, головные ударные волны и, так называемые, “магнитопаузы” (reconnection regions, bow shock and magnetopause). Одним из “удивительных”, если можно так выразиться, свойств магнитосферы Меркурия является наличие своеобразных вихрей, которые образуются благодаря взаимодействию внутреннего поля с солнечным ветром.

Многими исследователями предпринимались попытки построения аналитических моделей как “полного” магнитного поля, включающего в себя поле жидкого ядра [Margout et al., 2007] и поле магнитных масс, находящихся в коре Меркурия, так и отдельных составляющих магнитной индукции [Langlais et al., 2020]. Подход, основанный на разложении физических полей в ряд по сферическим гармоникам [Cain et al., 1989], является продуктивным в случае, когда измерения охватывают планету в глобальном масштабе [Wicht, Heyner, 2013]. Если же измерения более разрозненны и доступны только для отдельных областей планеты, то сферические гармоники также могут быть полезны, но при наложении дополнительных условий на коэффициенты разложения и выполнении процедуры регуляризации [Uno et al., 2009].

Магнитное поле Меркурия некоторые авторы [Langlais et al., 2020] пытаются описать как поле дипольных, квадрупольных и октупольных моментов. Перед выполнением интерпретации сигнал, соответствующий компоненте магнитного поля, проходит через узкополосный фильтр по методу Capon. Таким образом, восстанавливаются отдельные компоненты спектра магнитного поля, в том числе – компоненты поля литосферного происхождения. Нередко при решении подобного рода задач применяется метод дискретных дипольных моментов [Langlais et al., 2004; Langlais, Purucker, 2007; Hood et al., 2018]. Этот метод (дискретных точечных источников поля, сокращенно ESD) остается востребованным уже на протяжении порядка 40 лет [Mayhew, 1979].

Так называемые локальные методы тоже применяются при разработке аналитических моделей полей планет [Thébault et al., 2006; Hood et al., 2018].

Отличительные черты магнитного поля Меркурия накладывают ряд ограничений на построение модели магнитного динамо этой планеты, по сравнению с известными моделями палеодинамо Марса [Cao et al., 2014; Christensen, 2006].

В наших предыдущих работах [Kolotov et al., 2022; Stepanova, 2023] мы строили аналитические модели магнитного поля Меркурия (так называемое “магнитное изображение” Меркурия) на основе методики интегральных уравнений. Под магнитным “образом” или магнитным “изображением” Меркурия мы понимаем распределение магнитных диполей в коре планеты, эквивалентное по внешнему полю.

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

Миссия Messenger позволила также продвинуться в создании математических моделей гравитационного поля ближайшей к Солнцу планеты. Стоит отметить работу [Konopliv, 2020], в которой на основе данных станции Messenger были предложены три подхода (описанные как набор ограничений на возможные решения) к аппроксимации гравитационного поля Меркурия. Аналитические модели, описанные в указанной статье, согласованы с рассчитанной траекторией движения орбитальной станции: авторами предложен эффективный итерационный алгоритм уточнения модели гравитационного поля в зависимости от значений координат станции в планетоцентрической системе. Степень разрешения гравитационных данных варьировалась в широких пределах: от n = 12 в районе Южного полюса Меркурия до n = 154 в достаточно узкой области вблизи Северного полюса (такая детализация гравитационного поля стала возможной благодаря экстремально низким орбитам Messenger в конце периода работы станции). Через n здесь обозначен порядок гармоники при разложении поля в ряд по сферическим функциям.

В начале 2000-х был предложен весьма конструктивный подход, основанный на методе линейных интегральных представлений [Strakhov, Stepanova, 2002a; 2002b], который является обобщением в некотором смысле метода интегральных уравнений. С помощью различных версий метода линейных интегральных представлений были получены аналитические аппроксимации различных физических полей Земли и Марса [Gudkova et al., 2020; 2021; Salnikov et al., 2021; Раевский, Степанова, 2015а; 2015б; Степанова, 2020]. Метод линейных интегральных представлений позволяет находить устойчивые приближенные решения обратных задач геофизики даже при ограниченной информации об изучаемых элементах глобального или аномального поля, при наличии во входных данных помех различной природы. Основная идея метода заключается в представлении поля в виде простого и двойного слоев, распределенных на некоторых носителях (размерность носителей может принимать значения от единицы до четырех), а также в виде интегралов Радона и Фурье.

При совместной интерпретации данных по гравитационному и магнитному полям планет целесообразным представляется применение следующего подхода. Во-первых, данные спутникового зондирования должны быть отсортированы таким образом, чтобы распределение эквивалентных простого и двойного слоев можно было определить с минимально возможной ошибкой. Как этого добиться? Значения физических величин, измеренные космическими аппаратами в удаленных относительно поверхности планеты точках эллиптической орбиты, должны либо исключаться из выборки, по которой строится аналитическая модель поля, либо учитываться с весовыми множителями. Во-вторых, значения гравитационного и магнитного полей нужно пересчитать на оптимальной в некотором смысле сети наблюдений. Матрицы СЛАУ для компонент различных полей могут иметь как аналогичную, так и различную структуру: параметрический подход метода линейных интегральных представлений позволяет с высокой точностью восстанавливать пространственное распределение самых разнообразных сигналов, а также находить аналитическое продолжение гравитационного и магнитного полей вниз, в сторону источников, как было показано в наших предыдущих работах. Наша методика одновременной интерпретации “сырых” магнитных и синтетических гравитационных данных фактически является одним из способов цифровой обработки сигналов. Мы ищем такие распределения магнитных и гравитирующих масс в коре Меркурия, которые обеспечивают минимальное среднеквадратическое отклонение измеренных значений полей от их интегральных представлений, т.е. мы сравниваем искаженный интегральным представлением сигнал с его истинным изображением. В роли фильтра Винера выступает неизвестная плотность масс, в роли сигнала – значения производных разных порядков гравитационного потенциала, создаваемого точечными массами в заданных пунктах наблюдения. Результат свитки таких производных с неизвестной плотностью распределения – это и есть наблюдаемый сигнал. Возможны различные вариационные постановки, в которых учитываются корреляционные свойства сигнала и его изображения. Значения радиальных компонент гравитационного поля мы находили на трех сферических поверхностях, удаленных на 10, 50 и 100 км от поверхности планеты. При этом мы учитывали разложение гравитационного потенциала Меркурия по сферическим гармоникам до порядка 100. Магнитное поле было нами сначала проинтерпретировано с помощью метода интегральных уравнений, а затем пересчитано в точки той же самой сети, на которой рассматривалось гравитационное. В дальнейшем мы предполагаем построить оптимальную сеть наблюдений. Здесь имеется в виду такая схема расположения пунктов измерений гравитационного (или магнитного) поля, чтобы определитель матрицы СЛАУ, к которой редуцируется решение обратной линейной задачи, был отличен от нуля. Тогда решение СЛАУ находится однозначно, и этот факт позволяет существенно ограничить произвол в выборе эквивалентных по внешнему полю источников. Как правило, геофизические задачи отличаются тем, что СЛАУ плохо или очень плохо обусловленные – т.е. ранг матрицы СЛАУ меньше или существенно меньше числа неизвестных. Элементы матрицы СЛАУ в рамках метода S-аппроксимаций зависят от декартовых (локальный вариант) или сферических (региональный вариант) координат точек наблюдения. Если мы расположим точки наблюдения таким образом, чтобы определитель в нуль не обращался или, в крайнем случае, ранг матрицы был ненамного меньше числа неизвестных, то качество решения обратной линейной задачи по поиску распределения носителей масс на некоторой поверхности (в некотором объеме) значительно повысится.

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

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

$\left\{\begin{array}{l}\frac{\partial \overline{B}}{\partial t}=-\left\{\overline{v},\overline{B}\right\}+\eta \Delta \overline{B},\\ \text{div}B=\mathrm{0,}\end{array}\right.$

где: $\overline{B}(\overline{r},t)$ – это вектор магнитной индукции в некоторой области пространства М; $\overline{v}(\overline{r},t)$ – скорость жидкости; $\eta$ – малый безразмерный параметр, представляющий собой магнитную вязкость; $\left\{\overline{v},\overline{B}\right\}$ – это скобка Пуассона двух бездивергентных векторных полей. В трехмерном пространстве ее можно записать как $\left\{\overline{v},\overline{B}\right\}=-\text{rot}\left(\overline{v}\times \overline{B}\right).$ В более общей постановке, магнитное поле $\overline{B}(\overline{r},t)$ и поле скоростей $\overline{v}(\overline{r},t)$ определяются из системы уравнений

$\left\{\begin{array}{l}\frac{\partial \overline{B}}{\partial t}=-\left\{\overline{v},\overline{B}\right\}+\eta \Delta \overline{B},\\ \text{div}B=\mathrm{0,}\\ \frac{\partial \overline{v}}{\partial t}=-\left\{\overline{v},\overline{\nabla }\overline{v}\right\}+\text{rot}\overline{B}\times \overline{B}+\overline{v}\Delta \overline{v}-\nabla p.\end{array}\right.$

Здесь $\nabla p$ – градиент давления. Это – система уравнений магнитной гидродинамики.

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

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

Мы применили для совместной интерпретации гравитационных и магнитных данных Меркурия метод линейных интегральных представлений, основы которого изложены в работах [Strakhov, Stepanova, 2002a; 2002b]. В отличие от наших предыдущих работ, в данном исследовании рассматривались сразу два варианта вариационного подхода, а также выполнялось сравнение результатов моделирования гравитационного и магнитного полей в рамках вариационного и структурно-параметрического подходов. Последнее обстоятельство имеет принципиальное значение при цифровой обработке сигналов: необходимо подобрать в некотором смысле оптимальный фильтр, через который проходит искаженное “изображение” гравитационного или магнитного поля. Для поиска наилучшего фильтра необходимо формулировать условно-вариационные задачи таким образом, чтобы априорная информация об источниках поля (“истинном изображении”) учитывалась наиболее полно; кроме того, мелкие детали не должны отбраковываться при воздействии на искаженное поле. Метод линейных интегральных представлений и условно-вариационные постановки первого типа неоднократно рассматривались в наших работах [Страхов, Степанова, 2002; Степанова и др., 2020; Gudkova et al., 2020; 2021], мы хотели бы напомнить основные моменты при построении аналитических моделей физических полей.

Пусть задана совокупность величин $f_{i,\delta },$

$f_{i,\delta }=f_i+\delta f_i,$ $i=\mathrm{1,}\text{ 2, }\mathrm{...,}N,$ (1)

и при этом величины $f_i$ имеют представления вида

$f_i=\sum _{r=1}^R\underset{M_r}{\int }{\rho }_r\left(\xi \right)Q_r^{\left(i\right)}\left(\xi \right)d{\mu }_r\left(\xi \right),$ $i=\mathrm{1,}\text{ 2, }\mathrm{...,}N,$ (2)

где $Q_r^{\left(i\right)}\left(\xi \right),$ $i=\mathrm{1,}\text{ 2, }\mathrm{...,}N,$ $r=\mathrm{1,}\text{ 2, }\mathrm{...,}R,$ – заданные функции, удовлетворяющие условиям:

$\underset{M_r}{\overset{}{\int }}{\left(Q_r^{\left(i\right)}\left(\xi \right)\right)}^{2}d{\mu }_r\left(\xi \right)<+\infty ,$ (3)

и $M_r$ – области (конечные или бесконечные) или даже многообразия (поверхности, кривые) в R³, ${\mu }_r\left(\xi \right)$ – заданные на $M_r$ меры. Встречающиеся в подынтегральных выражениях в (2) функции ${\rho }_r\left(\xi \right)$ – неизвестные, подлежащие определению по заданным величинам $f_{i,\delta },$ $i=\mathrm{1,}\text{ 2, }\mathrm{...,}N,$ а также имеющейся априорной информации о векторе помехи $\delta f$ (N-векторе с компонентами $\delta f_i)$ и о свойствах искомых функций.

Основной смысл нахождения функций ${\rho }_r\left(\xi \right),$ т.е. построения интегральных представлений величин $f_i,$ см. (2), состоит в том, чтобы, получив приближенные (аппроксимирующие) представления:

${\rho }_r\left(\xi \right)\approx {\rho }_r^{\left(a\right)}\left(\xi \right)\text{, }\xi \in M_r,$ $r=\mathrm{1,}\text{ 2, }\mathrm{...,}R,$ (4)

далее использовать их для нахождения приближенных значений любых величин вида:

$p_s=\sum _{r=1}^R\underset{M_r}{\int }{\rho }_r\left(\xi \right)P_r^{\left(s\right)}\left(\xi \right)d{\mu }_r\left(\xi \right),$ s = 1, 2, ..., S,(5)

где $P_r^{\left(s\right)}\left(\xi \right)$ – заданные функции, по соотношениям:

$p_s\approx p_s^{\left(a\right)}=\sum _{r=1}^R\underset{M_r}{\int }{\rho }_r^{\left(a\right)}\left(\xi \right)P_r^{\left(s\right)}\left(\xi \right)d{\mu }_r\left(\xi \right),$ s = 1, 2, ..., S. (6)

Учет априорной информации о свойствах искомых функций ${\rho }_r\left(\xi \right)$ и вектора помехи $\delta f$ имеет важное значение при аналитическом моделировании любых физических полей.

Предположим, что существуют такие функции ${\rho }_r\left(\xi \right)>0:$

$\begin{array}{c}0\le m_{\text{min}}^{\left(r\right)}\le {\rho }_r\left(\xi \right)\le m_{\text{max}}^{\left(r\right)}<+\infty \text{,}\\ \xi \in M_r\text{,}\end{array}$ (7)

что выполняются неравенства:

$\underset{M_r}{\int }\frac{{\rho }_r^2\left(\xi \right)}{p_r\left(\xi \right)}d{\mu }_r\left(\xi \right)<+\infty ,$ $r=\mathrm{1,}\text{ 2, }\mathrm{...,}R.$ (8)

Неравенства (8) всегда имеют место, если одновременно выполняются два условия:

$\underset{\xi \in M_r}{\text{min}}{p}_r\left(\xi \right)>\mathrm{0,}$ $\underset{M_r}{\int }{\rho }_r^2\left(\xi \right)d{\mu }_r\left(\xi \right)<+\infty .$ (9)

Более специфическими свойствами функции ${\rho }_r\left(\xi \right)$ могут не обладать, что несущественно при нахождении интегральных представлений полей.

Но ситуация, когда $\underset{\xi \in M_r}{\text{min}}{\rho }_r\left(\xi \right)=0$ требует детального рассмотрения, поскольку большое практическое значение имеют случаи, когда $M_r$ – это конечные объемы или конечные куски поверхностей, на которых ставятся граничные условия:

${\left.\underset{}{{\rho }_r\left(\xi \right)}\right|}_{\partial M_r}=0.$ (10)

Эти условия удовлетворяются, если принять

${\left.\underset{}{p_r\left(\xi \right)}\right|}_{\partial M_r}\text{}=0.$ (11)

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

$0<{\delta }_{\text{min}}^2\le {||\delta f||}_{{\rm E}}^2\le {\delta }_{\text{max}}^2<+\infty$ (12)

и что априорно известным является равенство

$\left(f\text{, }\delta f\right)=0.$ (13)

В основной проблеме нахождения решения задачи построения аппроксимаций ${\rho }_r^{\left(a\right)}\left(\xi \right)$ функций ${\rho }_r\left(\xi \right)$ по заданному вектору $f_{\delta }$ (с компонентами $f_{i,\delta }),$ $f_{\delta }=f+\delta f,$ у которого компоненты f – полезного сигнала – имеют аналитические выражения (2), имеются два основных подхода: I) вариационный; II) структурно-параметрический.

Мы полагаем, что известна величина:

${\delta }^2={||\delta f||}_{{\rm E}}^2\text{\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\delta }_{\text{min}}^2\le {\delta }^2\le {\delta }_{\text{max}}^2.$ (14)

ВАРИАЦИОННЫЙ ПОДХОД

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

$\begin{array}{c}\Omega \left(\rho \right)=\sum _{r=1}^R\underset{M_r}{\int }\frac{{\rho }_r^2\left(\xi \right)}{p_r\left(\xi \right)}d{\mu }_r\left(\xi \right)=\underset{{\rho }_r^{}\left(\xi \right)}{\text{min}}\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},}\\ f_{i,\delta }-\sum _{r=1}^R\underset{M_r}{\int }{\rho }_r\left(\xi \right)Q_r^{\left(i\right)}\left(\xi \right)d{\mu }_r\left(\xi \right)=0\text{,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}i=\mathrm{1,}\mathrm{2,}\mathrm{...,}N\text{.}\end{array}$ (15)

Таким образом, в этом подходе минимизируется функционал $\Omega \left(\rho \right)$ при линейных – точечных – ограничениях.

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

$\begin{array}{c}\Omega \left(\rho \right)=\sum _{r=1}^R\underset{M_r}{\int }\frac{{\rho }_r^2\left(\xi \right)}{p_r\left(\xi \right)}d{\mu }_r\left(\xi \right)=\underset{{\rho }_r^{}\left(\xi \right)}{\text{min}}\text{,}\\ \sum _{i=1}^N{\left(f_{i,\delta }-\sum _{r=1}^R\underset{M_r}{\int }{\rho }_r\left(\xi \right)Q_r^{\left(i\right)}\left(\xi \right)d{\mu }_r\left(\xi \right)\right)}^2={\delta }^2\text{,}\end{array}$ (16)

где ${\delta }^2$ есть априорно заданная величина. Таким образом, во втором подходе минимизируется тот же самый функционал, но при одном квадратичном ограничении в виде равенства.

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

В случае безусловной вариационной задачи (15) (первый вариант вариационного метода) семейство безусловных вариационных задач таково:

$\begin{array}{c}\sum _{r=1}^R\underset{M_r}{\int }\frac{{\rho }_r^2\left(\xi \right)}{p_r\left(\xi \right)}d{\mu }_r\left(\xi \right)+\\ +2\sum _{i=1}^N{\lambda }_i\text{}\left(f_{i,\delta }-\sum _{r=1}^R\underset{M_r}{\int }{\rho }_r\left(\xi \right)Q_r^{\left(i\right)}\left(\xi \right)d{\mu }_r\left(\xi \right)\right)=\underset{{\rho }_r^{}\left(\xi \right)}{\text{min}}.\end{array}$ (17)

В (17) через ${\lambda }_i$ обозначены множители Лагранжа, с помощью которых учитываются фигурирующие в (15) линейные условия; в дальнейшем через $\lambda$обозначается вектор этих множителей

$\lambda ={\left({\lambda }_{\text{1}}\text{,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\lambda }_{\text{2}}\text{,\hspace{0.17em}\hspace{0.17em}}\mathrm{...,}{\lambda }_N\right)}^{\text{T}}.$ (18)

В случае же условной вариационной задачи (16) (второй вариант вариационного метода) семейство безусловных вариационных задач имеет вид:

$\begin{array}{c}\sum _{r=1}^R\underset{M_r}{\int }\frac{{\rho }_r^2\left(\xi \right)}{p_r\left(\xi \right)}d{\mu }_r\left(\xi \right)+\\ +\mu \sum _{i=1}^N{\left(f_{i,\delta }\text{}-\sum _{r=1}^R\underset{M_r}{\int }{\rho }_r\left(\xi \right)Q_r^{\left(i\right)}\left(\xi \right)d{\mu }_r\left(\xi \right)\right)}^2\text{}=\underset{{\rho }_r^{}\left(\xi \right)}{\text{min}},\end{array}$ (19)

где $\mu >0$ – множитель Лагранжа.

Методы вариационного исчисления позволяют сделать вывод, что в случае семейства безусловных вариационных задач (19) имеют место соотношения (здесь ${\tilde{\rho }}_r\left(\xi \text{;\hspace{0.33em}}\lambda \right)$ – решение задачи (19) при фиксированном векторе $\lambda ):$

${\tilde{\rho }}_r^{}\left(\xi \text{;}\lambda \right)=p_r\left(\xi \right)\sum _{i=1}^N{\lambda }_i{Q}_r^{\left(i\right)}\left(\xi \right),$ r = 1, 2, ..., R.(20)

Подчеркнем, что множители ${\lambda }_i$ в представлении функций ${\tilde{\rho }}_r\left(\xi \text{;\hspace{0.17em}}\lambda \right)$ одни и те же, т.е. не зависят от r.

В случае же семейства безусловных экстремальных задач (16) применение классических рецептов вариационного исчисления приводит к таким представлениям (здесь ${\tilde{\rho }}_r\left(\xi \text{;}\mu \right)$ – решение задачи (16) при фиксированном значении параметра $\mu ):$

${\tilde{\rho }}_r^{}\left(\xi \text{;}\mu \right)=\mu p_r\left(\xi \right)\sum _{i=1}^N\Delta f_{i,\mu }{Q}_r^{\left(i\right)}\left(\xi \right),$ r = 1, 2, ..., R, (21)

где через $\Delta f_{i,\mu }$ обозначены компоненты вектора

$\Delta f_{\mu }=f_{\delta }-{\tilde{f}}_{\mu },$ (22)

и

${\tilde{f}}_{i,\mu }=\sum _{r=1}^R\underset{M_r}{\int }{\tilde{\rho }}_r\left(\xi \text{;\hspace{0.17em}}\mu \right)Q_r^{\left(i\right)}\left(\xi \right)d{\mu }_r\left(\xi \right),$ i = 1, 2, ..., N, (23)

суть приближенные представления компонент $f_i$ вектора полезного сигнала. Отсюда ясно, что величины $\Delta f_{i,\mu }$ представляют собой оценки компонент вектора помехи $\delta f.$

Снова подчеркнем тот факт, что функции ${\tilde{\rho }}_r\left(\xi \text{;\hspace{0.33em}}\mu \right)$ зависят от одного и того же параметра $\mu ,$ который от r не зависит.

Далее. Ясно, что если использовать представления (20), то для нахождения вектора $\lambda$ множителей Лагранжа ${\lambda }_i$ получаем систему линейных алгебраических уравнений

$\text{A}\lambda =f_{\delta },$ (24)

в которой  является вектором входных данных задачи, см. (1), а A есть матрица с элементами

$\begin{array}{c}{a}_{ij}=\sum _{r=1}^R\underset{M_r}{\int }{p}_r^{}\left(\xi \right)Q_r^{\left(i\right)}\left(\xi \right)Q_r^{\left(j\right)}\left(\xi \right)d{\mu }_r\left(\xi \right)\text{,}\\ 1\le i\le N\text{,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}1\le j\le N\text{,\hspace{0.17em}\hspace{0.17em}}\end{array}$ (25)

т.е. обладающая свойством

$\text{A}={\text{A}}^{\text{T}}\ge 0.$ (26)

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

${||f_{\delta }-\text{A}{\tilde{\lambda }}_L||}_{{\rm E}}^2={\delta }^2,$ ${\delta }_{\min }^2\le {\delta }^2\le {\delta }_{\max }^2.$ (27)

Аналогично, если исходить из представлений (20)–(22), то для нахождения вектора $\Delta f_{\mu }$ получается система линейных алгебраических уравнений (при заданном значении параметра $\mu )$

$(\text{E}+\mu \text{A})\Delta f_{\mu }=f_{\delta },$ (28)

в которой A есть та же самая (N × N)-матрица, что и в системе (24), т.е. обладающая свойством (25) и имеющая элементы $a_{ij}$ по (26).

Из вышеизложенного можно сделать вывод о том, что и в рассматриваемом втором варианте вариационного метода нахождения приближенных представлений при наличии информации (13) имеем уравнение для нахождения значения параметра $\mu :$

${||\Delta f_{\mu }||}_{{\rm E}}^2={\delta }^2.$ ${\delta }_{\min }^2\le {\delta }^2\le {\delta }_{\max }^2.$ (29)

Система линейных алгебраических уравнений (28) при условиях (29) может применяться для нахождении фильтра Винера при решении задач цифровой обработки изображений. В принципе, с этой же целью можно рассматривать и постановки в рамках первого варианта вариационного подхода (15). Роль фильтра (или передаточной функции) играют функции (21) или (20), которые, как это будет показано ниже, могут представлять собой плотности простого и двойного слоев, распределенных на некоторых поверхностях (в частности, сферах). Изображениями считаются значения компонент магнитного или гравитационного поля в точках некоторой сети наблюдений (вопрос построения оптимальной сети наблюдений мы планируем рассмотреть в ближайшем будущем). А исходный сигнал – это значения различных производных гравитационного или магнитного потенциала в точках среды, специфическим образом связанных с точками сети наблюдений. Фактически, пункты наблюдений – это “зеркальные” изображения точечных источников (масс или диполей), расположенных под землей. “Зеркало” может быть как плоским (в случае локальных S-аппроксимаций), так и сферическим (региональные аппроксимации).

Далее заметим, что из уравнения (28) следует:

${\tilde{f}}_{\mu }=\mu \text{A}\Delta f_{\mu }.$ (30)

Иначе говоря, вектор оценки полезного сигнала ${\tilde{f}}_{\mu }$ выражается через вектор оценки погрешности (помехи) $\delta f.$ Но векторы f и $\delta f,$ вообще говоря, никак функционально не связаны. Поэтому дополнительным фактором, который может повысить устойчивость решения обратных задач к помехам во входных данных, мы считаем априорное задание функциональной зависимости векторов помехи и полезного сигнала, даже если зависимость подобного рода сводится к ортогональности векторов помехи и компоненты поля.

СТРУКТУРНО-ПАРАМЕТРИЧЕСКИЙ ПОДХОД

Переходим к рассмотрению второго, структурно-параметрического, подхода к нахождению приближенных представлений искомых функций ${\rho }_r^{}\left(\xi \right)$ по вектору экспериментальных данных о поле $f_{\delta }.$ Структурно-параметрический подход представляет собой чисто формальное обобщение основных приближенных представлений искомых функций ${\rho }_r^{}\left(\xi \right),$ возникающих в рамках вариационного подхода. Точнее говоря, соотношение (20) просто обобщается следующим образом:

$\begin{array}{c}{\tilde{\rho }}_r^{}\left(\xi \text{;}\stackrel{\circ }{{\lambda }_{}^{\left(r\right)}}\text{}\right)=p_r\left(\xi \right)\sum _{r=1}^N\stackrel{\circ }{{\lambda }_i^{\left(r\right)}}{Q}_r^{\left(i\right)}\left(\xi \right)\text{,}\\ r=\text{1,\hspace{0.33em}2, }\mathrm{...,}R\text{,\hspace{0.17em}\hspace{0.17em}}\end{array}$ (31)

(первый вариант структурно-параметрического подхода), а соотношение (21), соответственно, принимает вид:

$\begin{array}{c}{\tilde{\rho }}_r\left(\xi \text{;}{\mu }_r\right)={\mu }_r{p}_r\left(\xi \right)\sum _{i=1}^N\Delta f_{i,\mu }{Q}_r^{\left(i\right)}\left(\xi \right)\text{,}\\ r=\text{1,\hspace{0.33em}2, }\mathrm{...,}R\text{,\hspace{0.17em}\hspace{0.17em}}\end{array}$ (32)

где теперь

$\mu ={\left({\mu }_{\text{1}}\text{,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\mu }_{\text{2}}\text{, }\mathrm{...,}{\mu }_R\right)}^{\text{T}},$ (33)

и вектор $\Delta f_{\mu }$ (с компонентами $\Delta f_{i,\mu })$ зависит уже не от скаляра $\mu ,$ а от вектора $\mu$ с R-компонентами, определенными по (33).

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

Нетрудно показать, что в случае первого варианта структурно-параметрического метода получаем систему линейных алгебраических уравнений:

$\stackrel{\circ }{\text{A}}\stackrel{\circ }{\lambda }=f_{\delta },$ (34)

в которой $\stackrel{\circ }{\lambda }$ есть блочный вектор вида:

$\stackrel{\circ }{\lambda }=\begin{array}{ccc}& \stackrel{\circ }{{\lambda }_{}^{\left(1\right)}}& \\ & \stackrel{\circ }{{\lambda }_{}^{\left(2\right)}}& \\ & & \\ & \stackrel{\circ }{{\lambda }_{}^{\left(R\right)}}& \end{array}\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},\hspace{1em}\hspace{1em}\hspace{0.17em}}\stackrel{\circ }{{\lambda }_{}^{\left(r\right)}}=\begin{array}{ccc}& \stackrel{\circ }{{\lambda }_1^{\left(r\right)}}& \\ & \stackrel{\circ }{{\lambda }_2^{\left(r\right)}}& \\ & & \\ & \stackrel{\circ }{{\lambda }_N^{\left(r\right)}}& \end{array}\text{\hspace{0.17em}\hspace{0.17em},}$ (35)

а матрица $\stackrel{\circ }{\text{A}}$ суть блочная матрица следующего вида:

$\stackrel{\circ }{\text{A}}=\begin{array}{cccc}\stackrel{\circ }{{\text{A}}^{\text{(1)}}}& \stackrel{\circ }{{\text{A}}^{\text{(2)}}}& \begin{array}{cc}& \end{array}& \stackrel{\circ }{{\text{\hspace{0.33em}A}}^{\text{(}R\text{)}}}\end{array},$ (36)

в которой блоки $\stackrel{\circ }{{\text{A}}^{\text{(}r\text{)}}}$ имеют элементы

$\begin{array}{c}\stackrel{\circ }{a_{ij}^{\left(r\right)}}\text{}=\underset{M_r}{\int }{p}_r\left(\xi \right)Q_r^{\left(i\right)}\left(\xi \right)Q_r^{\left(j\right)}\left(\xi \right)d{\mu }_r\left(\xi \right)\text{,}\\ r=\text{1, 2, }\mathrm{...,}R\end{array}$ (37)

и обладают свойством:

$\begin{array}{c}\stackrel{\circ }{{\text{A}}^{\text{(}r\text{)}}}\text{}=\stackrel{°}{{\text{A}}^{\text{(}r\text{),T}}}\ge 0\text{\hspace{0.17em},}\\ r=\text{1, 2, }\mathrm{...,}R.\end{array}$ (38)

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

Структурно-параметрический подход, примененный к решению обратных задач геофизики в рамках метода линейных интегральных представлений, обладает рядом замечательных свойств: он позволяет находить согласованные с априорной информацией о сигнале распределения эквивалентных источников гравитационного и магнитного полей при рассмотрении вариационных постановок в пространстве функций, интегрируемых с квадратом, а не в пространстве функций с равномерной нормой (или дифференцируемых функций с ограничениями нормы самой функции и ее производных). Это – принципиальный момент при выполнении исследований в области линейных и (особенно!) нелинейных обратных задач геофизики и при восстановлении источников сигналов любой природы по искаженным, в общем случае, изображениям. Зависимость вектора параметров в (37) от структурного элемента (выражающаяся в наличии индекса “r” у компонент вектора решения) позволяет построить простые аналитические модели различных физических полей. Мы можем по измеренным на одной и той же сети точек значениям различных компонент гравитационного поля, магнитного поля и другим данным об объекте исследования воссоздать источники этих полей с приемлемой точностью в случае наличия априорной информации об источниках, а также предложить метрологические аппроксимации полей, если никакой количественной информации о мощности и геометрии гравитирующих и магнитных масс нет (а именно так и бывает при изучении планет).

Если известны компоненты магнитного или гравитационного поля (например, первая производная потенциала по r на некотором рельефе), то можно представить компоненты поля в виде суммы простого и двойного слоев, создаваемых совокупностью горизонтальных плоскостей, расположенных ниже заданного рельефа.

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

Представим себе идеализированный Меркурий как внутренность сферы радиуса $R_0.$ Реальный Меркурий тогда может трактоваться как область трехмерного пространства, ограниченная замкнутой кусочно-непрерывной поверхностью S, незначительно отклоняющейся от сферы радиуса $R_0.$Предполагается, что приближенные значения функции $G\left(x\right),$ гармонической вне этой сферы, заданы в произвольном множестве точек $x^{\left(i\right)}\text{},i=\mathrm{1,}\mathrm{2,}\dots N,$ на поверхности S:

$\begin{array}{c}{f}_{i,\delta }=f_i+\delta f_i,\\ f_i=G\left(x^{\left(i\right)}\right).\end{array}$ (39)

Так как G(x) гармонична при $r>R_0,$ она имеет следующее интегральное представление:

$\begin{array}{c}G\left(x\right)=\frac{R_0^2}{4\pi }\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\pi }{\int }}\frac{\sigma (\tilde{\vartheta },\tilde{\phi })\sin \tilde{\vartheta }d\tilde{\vartheta }d\tilde{\phi }}{R(\xi -x)}+\\ +\frac{R_0^2}{4\pi }\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\pi }{\int }}\frac{w(\tilde{\vartheta },\tilde{\phi })(R_0-r\cos {\vartheta }^{\text{'}})\sin \tilde{\vartheta }d\tilde{\vartheta }d\tilde{\phi }}{R^3(\xi -x)},\\ R(\xi -x)={\left(R_0^2-2R_{0}r\cos {\vartheta }^{\text{'}}+r^2\right)}^{1/2},\\ \xi =(R_0\cos \tilde{\phi }\sin \tilde{\vartheta },R_0\sin \tilde{\phi }\sin \tilde{\vartheta },R_0\cos \tilde{\vartheta }),\\ x=(r\cos \phi \sin \vartheta ,r\sin \phi \sin \vartheta ,r\cos \vartheta ),\\ \cos {\vartheta }^{\text{'}}=\sin \vartheta \sin \tilde{\vartheta }\cos (\phi -\tilde{\phi })+\cos \vartheta \cos \tilde{\vartheta }.\end{array}$ (40)

Функция $\sigma (\tilde{\vartheta },\tilde{\phi })$ в (40) называется плотностью простого слоя, распределенного на сфере радиуса $R_0,$ $w(\tilde{\vartheta },\tilde{\phi })$ – плотность двойного слоя, распределенного на той же сфере; $R(\xi -x)$ – расстояние между текущей точкой сферы $\xi$ и точкой наблюдения x. Дифференцирование правой части (40) по различным координатам точки наблюдения G(x), даже если эти производные (например, $\frac{\partial V}{\partial r}$ – радиальная производная гравитационного потенциала) не являются гармоническими функциями.

Формула (40) является одной из форм интегрального представления гармонической функции вне сферы радиуса R [Владимиров, 1981]:

$V\left(x\right)=\frac{1}{4\pi }\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\pi }{\int }}\left(\frac{1}{R(\xi -x)}\frac{\partial V\left(\xi \right)}{\partial \overrightarrow{n}}-V\left(\xi \right)\frac{\partial }{\partial \overrightarrow{n}}\left(\frac{1}{R(\xi -x}\right)\right)\sin \tilde{\vartheta }d\tilde{\vartheta }d\tilde{\phi },$ (41)

где $\overrightarrow{n}$ – внешняя нормаль к единичной сфере (мы можем написать $\frac{\partial }{\partial r},$ r – радиальная координата радиус-вектора x).

Условно-вариационная задача для $\sigma$ и w [Strakhov, Stepanova, 2002b] в рамках второго варианта вариационного подхода (см. (17)) имеет следующее решение:

$\begin{array}{c}\sigma \left(\hat{\xi }\right)=\sum _{i=1}^N{\lambda }_i{Q}_i^{\left(1\right)}\left(\xi \right),\\ w\left(\hat{\xi }\right)=\sum _{i=1}^N{\lambda }_i{Q}_i^{\left(2\right)}\left(\xi \right),\hat{\xi }=(\tilde{\vartheta },\tilde{\phi }),\xi =(R,\tilde{\vartheta },\tilde{\phi }),\end{array}$ (42)

где $Q_i^{\left(1\right)}\left(\xi \right)$ и $Q_i^{\left(2\right)}\left(\xi \right)$ представляются в виде

$Q_i^{\left(1\right)}\left(\xi \right)=\frac{1}{R(\xi -x^{\left(i\right)})},Q_i^{\left(2\right)}\left(\xi \right)=\frac{R_0-r_i\cos {{\vartheta }^{\text{'}}}_i}{R^3(\xi -x^{\left(i\right)})}.$ (43)

Величины ${\lambda }_i$ являются компонентами N-вектора $\lambda ,$ дающего решение системы линейных уравнений:

$A\lambda =f_{\delta },f_{\delta }=f+\delta f,$ (44)

где $f_{\delta }$ имеет компоненты $f_{i,\delta }$ (39) и элементы матрицы $A=A^T\ge 0$ записываются как

$\begin{array}{c}{a}_{ij}=\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\pi }{\int }}\left[Q_i^{\left(1\right)}\left(\xi \right)Q_j^{\left(1\right)}\left(\xi \right)+Q_i^{\left(2\right)}\left(\xi \right)Q_j^{\left(2\right)}\left(\xi \right)\right]\sin \tilde{\vartheta }d\tilde{\vartheta }d\tilde{\phi }=\\ =\frac{2\pi }{{\left(h_i{h}_j\right)}^{1/2}{r}_i{r}_j}\left(1+\frac{0.25}{r_i{r}_j{h}_i{h}_j}\right)F\left(2\text{arctg}\left(\sqrt{h_i{h}_j}\right),\sqrt{0.5\cdot \left(1\text{}+\cos {{\alpha }_i}_j\right)}\right)-\\ -\frac{\pi \left(3{\left(h_i{h}_j\right)}^2-4h_i{h}_j\cos {{\alpha }_i}_j+1\right)}{{r_i}^2{r}_j^2{h}_i{h}_j{\left(\sqrt{\left.1-2h_i{h}_j\cos {{\alpha }_i}_j\right)+{\left(h_i{h}_j\right)}^2}\right)}^3}.\end{array}$ (45)

Здесь $h_i=\frac{R_0}{r_i},h_j=\frac{R_0}{r_j};$ ${\alpha }_{ij}$ – угол между векторами $x_i$ и $x_j.$ Мы предполагаем, что вектор $x_i$ параллелен оси Oz , а вектор $x_j$ лежит в плоскости zOx (мы всегда так можем выбрать систему координат). Тогда указанные вектора будут иметь координаты:

$x_i=(r_i,\mathrm{0,}0)$ и $x_j=(r_j\sin {\alpha }_{ij},\mathrm{0,}{r}_j\cos {\alpha }_{ij}).$

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

Если ограничиться представлением элементов магнитного поля Меркурия в виде потенциала простого слоя, то элементы матрицы системы (45) принимают тогда вид:

$a_{ij}=\frac{2\pi }{{\left(h_i{h}_j\right)}^{1/2}{r}_i{r}_j}F\left(2\text{arctg}\left(\sqrt{h_i{h}_j}\right),\sqrt{0.5\cdot \left(1+\cos {{\alpha }_i}_j\right)}\right).$ (46)

Если применяется первый вариант вариационного подхода (см. (16)), то матрица системы (45) характеризуется элементами:

$\begin{array}{c}{a}_{ij}={\delta }_{ij}+\mu \cdot \frac{2\pi }{{\left(h_i{h}_j\right)}^{1/2}{r}_i{r}_j}\times \\ \times F\left(2\text{arctg}\left(\sqrt{h_i{h}_j}\right),\sqrt{0.5\cdot \left(1+\cos {{\alpha }_i}_j\right)}\right).\end{array}$ (47)

Значения параметра  можно вычислять по рекуррентной формуле:

${\mu }^{(k+1)}=\alpha {\mu }^{\left(k\right)},k=\mathrm{0,}\mathrm{...,}\infty ;{\mu }^{\left(0\right)}=0.5;0<\alpha <1.$ (48)

Однако параметры метода можно определять и иначе: можно считать, что $\mu$ – это еще одна неизвестная величина, которая восстанавливается из решения некоторой системы алгебраических уравнений (подчеркнем, соответствующая система уравнений является нелинейной).

РЕЗУЛЬТАТЫ МАТЕМАТИЧЕСКОГО ЭКСПЕРИМЕНТА

Данные миссии Messenger за несколько дней полета были интерпретированы с помощью трех различных модификаций метода линейных интегральных представлений. Точек в наборе было 10 000. В файлах, содержащих “сырые данные”, указывались декартовы координаты точек наблюдения в километрах, при этом начало системы координат совпадает с центром масс Меркурия. Носители простого и двойного слоев при аппроксимации лапласовых полей располагались в коре Меркурия, т.е. на расстоянии от 0.1 до 100 км от поверхности планеты.

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

Структурно-параметрический подход подразумевает следующее: для каждого из носителей определяется свой вектор решения [Salnikov et al., 2021]. Во всех случаях системы линейных алгебраических уравнений (СЛАУ), к которым редуцируется решение обратной задачи по восстановлению магнитного поля Меркурия, решаются с помощью метода регуляризации Холецкого (БМХР) и усовершенствованного блочного метода решения СЛАУ (УБМ) [Раевский, Степанова, 2015а; 2015б; Gudkova et al., 2020; Степанова и др., 2020; Salnikov et al., 2021]. При этом мы полагали, что Меркурий представляет собой шар, радиуса $R_0=2439.0$ км. Результаты аппроксимации представлены в таблице.

 

Модифицированные S-аппроксимации z-компоненты магнитного и гравитационного полей Меркурия по данным Messenger

№	$\frac{B_r}{N},\frac{g}{N}$	$R_1,R_2,$ км R	Метод решения СЛАУ	${\sigma }_{\min },$ нТл или мГал	${\sigma }_{\max },$ нТл или мГал	${\sigma }_0,$ нТл или мГал	$\frac{\Delta }{t}$
1	$\begin{array}{c}{B}_r\\ 10000\end{array}$	2400, 2380	БМХР	0.012	0.024	0.015	$\begin{array}{c}1.1\times {10}^{-3}\\ 34:15\end{array}$
2	$\begin{array}{c}{B}_r\\ 7000\end{array}$	2400, 2380, 2489	БМХР	0.0027	0.0039	0.0029	$\begin{array}{c}1.4\times {10}^{-6}\\ 15:34\end{array}$
3	$\begin{array}{c}{B}_r\\ 7000\end{array}$	2400, 2380, 2539	УБМ	0.001	0.0015	0.0014	$\begin{array}{c}1.5\times {10}^{-8}\\ 11:12\end{array}$
4	$\begin{array}{c}g\\ 7000\end{array}$	2400, 2380, 2449	УБМ	0.001	0.0015	0.0011	$\begin{array}{c}1.24\times {10}^{-8}\\ 45:22\end{array}$
5	$\begin{array}{c}g\\ 7000\end{array}$	2400, 2380, 2489	УБМ	0.01	0.015	0.014	$\begin{array}{c}1.22\times {10}^{-4}\\ 20:36\end{array}$
6	$\begin{array}{c}g\\ 7000\end{array}$	2400, 2380, 2539	УБМ	0.001	0.0015	0.0011	$\begin{array}{c}1.16\times {10}^{-7}\\ 31:15\end{array}$

Примечания: $\sigma =\frac{{||Ax-f_{\delta }||}_E}{\sqrt{N}}$ – среднеквадратическое отклонение; ${\sigma }_{\max }=\sqrt{\frac{{\delta }_{\max }^2}{N}},$ ${\sigma }_{\min }=\sqrt{\frac{{\delta }_{\min }^2}{N}};$ $\Delta =\frac{{||Ax-f_{\delta }||}_E}{{||f_{\delta }||}_E}$ – показатель качества решения; ${\sigma }_0$ – среднеквадратическое отклонение, полученное в результате решения СЛАУ; t – время в минутах и секундах.

 

Орбита космической миссии в некоторые моменты времени достаточно далеко уходила от поверхности Меркурия (на расстояние до 0.6 среднего радиуса планеты), поэтому для выделения из “сырых” данных составляющих магнитной индукции, генерируемых токами в жидком ядре и коре (так называемого внутреннего магнитного поля Меркурия), можно воспользоваться приближением “the thin shell approximation” [Toepfer et al., 2021]. При таком подходе точки наблюдения должны находиться в пределах “тонкой” (по сравнению с некоторыми параметрами, характеризующими топологию планеты) оболочки, окружающей Меркурий. Полоидальное и тороидальное магнитные поля, создаваемые токами в плазме вокруг Меркурия, “исчезают” в этом случае. Мы проводили математический эксперимент, не выходя за рамки модели “тонкой оболочки”: каждый отрезок из некоторого набора данных, полученных при движении спутника, находится в пределах сферической оболочки, толщина которой составляет приблизительно 0.1 радиуса Меркурия, т.е. 240 км. Процесс создания математической модели магнитного поля по “сырым” спутниковым данным подразумевает прохождение следующих этапов.

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

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

Далее, в точках той же самой сети наблюдений, для которой были вычислены значения радиальной компоненты магнитного поля по известной формуле для разложения гравитационного потенциала в ряд по полностью нормированным присоединенным функциям Лежандра [Gudkova et al., 2020], нами были выполнены расчеты радиальной компоненты гравитационного поля Меркурия.

Данные гравитационного поля и топографии Меркурия, представленные в виде разложения относительно центра масс по нормализованным сферическим функциям, можно найти на сайте Системы Планетных данных (http://pds-geosciences.wustl.edu):

$V\left(r,\phi ,\lambda \right)=\frac{GM}{r}\left[1+\text{}\sum _{n=2}^{ng}\sum _{m=0}^n{\left(\frac{R_e}{r}\right)}^n\text{}\left(C_{gnm}\text{}\cos m\lambda +S_{gnm}\text{}\sin m\lambda \right)P_{nm}\text{}\left(\sin \phi \right)\right],$ (49)

$\begin{array}{c}{R}_{рельеф}(r,\phi ,\lambda )=\\ =R_0+\sum _{n=1}^{nt}\sum _{m=0}^n(C_{tnm}\cos m\lambda +S_{tnm}\sin m\lambda )P_{nm}\left(\sin \phi \right),\end{array}$ (50)

где: V(r, ϕ, λ) – гравитационный потенциал; r – расстояние от точки наблюдения до начала координат; ϕ и λ – ареоцентрические широта и долгота; М – масса планеты; G – гравитационная постоянная; R_e – экваториальный радиус (нормирующий радиус); R₀ – средний радиус планеты; C_gnm и S_gnm – нормализованные гравитационные моменты; C_tnm и S_tnm – гармонические коэффициенты в разложении рельефа (в метрах); ng и nt – степени гармоник до которых проведено разложение гравитационного поля и топографии по сферическим функциям соответственно; P_nm – нормализованные функции Лежандра степени n и азимутального числа m, удовлетворяющие соотношению:

$\underset{0}{\overset{2\pi }{\int }}\underset{-1}{\overset{1}{\int }}{P}_{nm}^2\left(\sin \phi \right)\left\{\begin{array}{c}{\cos }^2\left(m\lambda \right)\\ {\sin }^2\left(m\lambda \right)\end{array}\right.d\left(\sin \phi \right)d\lambda =4\pi .$

При дифференцировании гравитационного потенциала по радиусу получается выражение для радиальной компоненты гравитационного поля. Мы рассматривали разложения гравитационного поля Меркурия до 100 гармоники включительно. С помощью двух вариантов вариационного подхода и в рамках структурно-параметрического подхода нами определялись эквивалентные по внешнему гравитационному полю распределения масс, залегающие на нескольких сферических поверхностях в коре Меркурия. Для “чистоты” эксперимента параметры сфер выбирались такими же, как и в случае “сырых” данных по магнитному полю. Необходимо подчеркнуть следующее: метод линейных интегральных представлений в различных вариантах и при применении различных подходов позволяет единообразно описать процесс построения математической модели любого сигнала. В настоящей работе мы исследуем возможность решения обратных линейных задач для стационарных сигналов, но, как указывалось нами в работе [Степанова и др., 2022], эффективные алгоритмы интерпретации данных можно строить и для дифференциальных операторов, зависящих от времени, т.е. для нестационарных полей. Ключевым вопросом при одновременной интерпретации разнородных данных в рамках описанного в настоящей работе подхода является следующий: как влияют параметры сети наблюдений и скорость убывания сигнала в зависимости от расстояния от источника на точноcть аппроксимации компонент поля и на точность решения линейных обратных задач? Что понимается под точностью решения линейной обратной задачи при неизвестных распределениях масс (как магнитных, так и гравитационных)? Мы предлагаем следующий критерий условной точности решения: если по разным выборкам аппроксимированного поля получаются мало отличающиеся распределения эквивалентных источников, то такое решение обратной линейной задачи можно назвать адекватным реальным данным. Количественно можно описать предложенный критерий качества решения так:

${\Delta }^{\text{'}}=\frac{{||x^{\left(k\right)}-x^{(k-1)}||}_E}{{||x^{(k-1)}||}_E}\le \epsilon ,k=\mathrm{0,}\mathrm{1,}\mathrm{....}$ (51)

В (51) под $x^{\left(k\right)}$ понимается плотность распределения потенциала простого или двойного слоя, полученная по k-й выборке значений поля. В наших математических экспериментах мы восстанавливали фиктивные источники по значениям гравитационного (синтезированного!) и магнитного (также синтезированного, но другим способом – с помощью модифицированных S-аппроксимаций) полей на сферах, отстоящих от поверхности Меркурия на 10, 50 и 100 км. Подчеркнем, что мы не “удалялись” от поверхности планеты на слишком большое расстояние (максимум – на 100 км).

В таблице приводятся результаты аппроксимации магнитного и гравитационного полей по пересчитанным на указанные сферы данным, а также результат решения СЛАУ по “сырым” магнитным данным (пример 1). Примеры 2–3 в таблице соответствуют аппроксимации магнитного поля по синтезированным данным на расстоянии 50 и 100 км от поверхности планеты в рамках структурно-параметрического подхода; примеры 4–6 демонстрируют реализацию структурно-параметрического подхода при нахождении модифицированных региональных S-аппроксимаций синтезированного на расстоянии 10, 50 и 100 км, соответственно, гравитационного поля.

На рис. 1a изображено магнитное поле (r-компонента) Меркурия по данным миссии Messenger; на рис. 1б представлено магнитное поле вблизи поверхности, построенное с помощью региональных S-аппроксимаций в рамках структурно-параметрического подхода (пример 3 в таблице. На рис. 2 представлена проекция орбит станции Messenger на плоскость экватора.

 

Рис.  1. (а) – Bz-компонента магнитного поля Меркурия по данным миссии Messenger; (б) – карта изолиний магнитного поля Меркурия вблизи поверхности, построенного с помощью региональных S-аппроксимаций.

 

Рис.  2. (а) – Схематическое изображение проекций траекторий пролета станции Messenger на плоскость экватора Меркурия; (б) – схематическое изображение траекторий пролета станции Messenger как поверхности в трехмерном пространстве.

 

На рис. 3 и рис. 4 изображены распределения эквивалентных магнитных масс, восстановленных по синтезированному магнитному полю на расстоянии 50 и 100 км от поверхности планеты. Синтез радиальной компоненты магнитного поля осуществлялся по полученным решениям обратной задачи для “сырых” данных, измеренных на эллиптической орбите. Мы рассчитали модельное гравитационное поле по формулам (49), (50), значения радиальной компоненты гравитационного поля Меркурия приведены на рис. 5. На рис. 6 показано распределение гравитационных масс, определенных по синтезированному гравитационному полю (см. формулы (49)–(51)) на сферической поверхности, отстоящей от поверхности Меркурия на 50 км соответственно. Из рис. 3, рис. 4 видно, что характер распределения эквивалентных носителей практически не меняется в зависимости от сферы, на которой вычислялось синтезированное поле. Характеристика “точности” эквивалентного решения (см. (51)) практически не менялась при расчетах на сферах, удаленных от поверхности планеты на расстояние не более 100 км. Для гравитационного поля она составила 1–5%, для магнитного – 3–7%.

 

Рис.  3. Распределение магнитных масс, найденное по аналитически продолженному на 50 км вверх магнитному полю Меркурия.

 

Рис.  4. Распределение магнитных масс, найденное по аналитически продолженному на 100 км вверх магнитному полю Меркурия.

 

Рис.  5. Модельное гравитационное поле Меркурия.

 

Рис.  6. Распределение эквивалентных гравитационных источников, определенное по полю на высоте 50 км от поверхности Меркурия.

 

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

Значение показателя качества решения может быть меньше, чем относительная погрешность данных наблюдений. Это происходит в тех случаях, когда СЛАУ относительно хорошо обусловлена и решение СЛАУ находится с помощью регуляризованных методов практически сразу, число итераций по параметру регуляризации не превышает 2–3. Однако при построении аппроксимаций по “сырым” данным, такое случается сравнительно редко, и поэтому алгоритм нахождения пробных решений [Strakhov, Stepanova, 2002a; 2002b], удовлетворяющих ограничениям по невязке как сверху, так и снизу, остается актуальным.

ЗАКЛЮЧЕНИЕ

Предлагаемая в настоящей работе методика одновременной интерпретации данных о магнитном и гравитационном полях может быть полезна при цифровой обработке изображений, поскольку нахождение эквивалентных по внешнему полю плотностей простого и двойного слоев в рамках второго варианта вариационного подхода аналогично определению фильтра Винера при анализе сигнала и его изображения. Различие заключается в том, что измеренные вне планеты поля являются изображениями источников, расположенных внутри планеты, т.е. область определения изображения не совпадает с областью определения исходного сигнала. Но не совпадают только значения z-координат источника и изображения: интегрирование по горизонтальным координатам в локальном варианте метода S-аппроксимаций или по сферам в региональной версии того же метода выполняется одинаково для сигнала и его изображения. Роль фильтра в нашем случае играет плотность простого или двойного слоя, а роль источника сигнала – производная потенциала точечного источника гравитационного или магнитного поля, т.е. функция Грина соответствующей краевой задачи. Это становится очевидным, если внимательно присмотреться к интегральным представлениям и виду плотностей простого и двойного слоев (41)–(43).

Мы упомянули во Введении так называемый метод ESD [Emilia, 1973; Langlais et al., 2004]. Существуют модификации этого метода, учитывающие зависимость магнитных свойств среды от времени. В основе этого метода лежит представление полного поля Меркурия в виде суперпозиции полей, создаваемых глубоко залегающими магнитными диполями. Количество таких диполей может достигать нескольких сот миллионов. Авторы метода предполагают, что источники залегают глубоко в недрах планеты – вблизи ядра, регуляризирующие алгоритмы при применении метода, в основном, не рассматриваются. Помимо этого, авторы метода ESD полагают, что разложение магнитного поля планеты по сферическим гармоникам дает надежные результаты лишь при охвате всей планеты в целом, т.е. локальные и региональные неоднородности определить в рамках указанного метода не представляется возможным.

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

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

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

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

В настоящей работе мы построили аналитические аппроксимации Br-компоненты магнитного поля и радиальной производной гравитационного потенциала Меркурия в сферической системе координат, начало которой жестко связано с центром масс планеты. Эллиптичность планеты не учитывалась. Измерения магнитного поля выполнялись с помощью аппаратуры, установленной на станции Messenger, при ее пролетах вокруг Меркурия (мы постарались учесть результаты наблюдений на не слишком больших расстояниях от поверхности планеты). Также была предпринята попытка аналитического продолжения магнитного поля, измеренного на орбите спутником, вниз, в сторону источников. Таким образом, мы воссоздали поле на поверхности планеты и на расстояниях, приблизительно равных 10, 50 и 100 км от Меркурия. Относительная точность математического моделирования поля на спутниковой орбите оказалась достаточно высокой – она не превышала $1.1\times {10}^{-3}.$ Магнитометры, установленные на Messenger, могли измерять поле с точностью до 0.0004 нТл, но при оценке показателя качества решения, $\Delta =\frac{{||Ax-f_{\delta }||}_E}{{||f_{\delta }||}_E},$ необходимо учитывать множество факторов: ошибки в определении координат космического аппарата, систематические ошибки приборов, взаимодействие магнитного поля планеты и солнечного ветра и т.п. Реальная точность измерений, по оценкам разных авторов, не превышает 0.05 нТл. По нашему мнению, основным моментом при разработке методики математического моделирования физических полей Земли и планет должно быть соответствие полученных результатов реальным, “сырым”, данным. Аппроксимированное поле не должно “распадаться” при приближении к поверхности планеты. Ввиду существенной неоднозначности решения обратных задач геофизики, нахождение эквивалентных по внешнему полю распределений масс является одним из возможных способов интерпретации разнородных и разноточных данных.

ФИНАНСИРОВАНИЕ РАБОТЫ

Работа выполнена при поддержке РНФ (проект № 23-41-00002).

Авторлар туралы

I.  Stepanova

Schmidt Institute of Physics of the Erath, Russian Academy of Sciences

Хат алмасуға жауапты Автор.
Email: tet@ifz.ru
Ресей, Moscow

A.  Yagola

Moscow State University

Email: tet@ifz.ru
Ресей, Moscow

D.  Lukyanenko

Moscow State University

Email: tet@ifz.ru
Ресей, Moscow

I.  Kolotov

Moscow State University

Email: tet@ifz.ru
Ресей, Moscow

Әдебиет тізімі

Қосымша файлдар