1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 60, № 6, 2020

Журнал вычислительной математики и математической физики, 2020, T. 60, № 6, стр. 1027-1034

Восстановление магнитной восприимчивости с использованием полных магнито-градиентных данных

Я. Ван 1, 2, 3, *, И. И. Колотов 4, **, Д. В. Лукьяненко 4, ***, А. Г. Ягола 4, ****

1 Китайская академия наук, Институт геологии и геофизики, Ключевая лаборатория исследований нефтяных ресурсов
100029 Пекин, Китай

2 Университет Китайской академии наук
100049 Пекин, Китай

3 Китайская академия наук, Институт наук о Земле
100029 Пекин, Китай

4 МГУ, физический факультет, кафедра математики
119992 Москва, Ленинские горы, Россия

* E-mail: yfwang@mail.iggcas.ac.cn
** E-mail: kolotovigor@list.ru
*** E-mail: lukyanenko@physics.msu.ru
**** E-mail: yagola@physics.msu.ru

Поступила в редакцию 28.10.2019
После доработки 28.10.2019
Принята к публикации 11.02.2020

Полный текст (PDF)

Аннотация

В работе рассмотрены особенности решения обратной задачи восстановления магнитной восприимчивости с использованием полных тензорных магнитно-градиентных данных. Данная задача сводится к решению системы двух трехмерных интегральных уравнений Фредгольма I рода, одно из которых связывает магнитную восприимчивость ограниченного тела с индуцированным им магнитным полем, другое – с полным тензором градиентов компонент магнитной индукции. На серии модельных примеров продемонстрировано, что использование полного тензора градиентов компонент магнитной индукции существенным образом улучшает качество восстанавливаемой функции, определяющей магнитную восприимчивость. Библ. 23. Фиг. 1.

Ключевые слова: магнитостатика, магнитная восприимчивость, полный тензор градиентов компонент магнитной индукции, обратная задача.

ВВЕДЕНИЕ

Математические модели в магнитостатике, связывающие плотность магнитного момента ограниченного тела и его магнитную индукцию, достаточно хорошо известны [1]–[4]. Менее известны математические модели магнитостатики, связывающие объемную (или поверхностную) плотность магнитного момента тела с полным тензором градиентов компонент магнитной индукции (ГКМИ), измеряемым в некоторой удаленной от тела области [5]–[10]. Китайским участником работы вместе со своими коллегами из Института геологии и геофизики Китайской академии наук было установлено [11], [12], что тензор градиентов компонент магнитной индукции в таких моделях обладает большей чувствительностью к тонкой структуре распределения магнитного момента, чем сама магнитная индукция. Поэтому измерения тензора ГКМИ представляются более перспективными для интерпретации магнитных полей с помощью решения обратных задач.

Модели, включающие тензор ГКМИ, могут быть рассмотрены в различных вариантах: линейных и нелинейных, непрерывных и дискретных, для двумерных и трехмерных магнитных полей. Все эти модели порождают обратные задачи нахождения плотности магнитного момента по тензору ГКМИ в различных постановках задачи. Одной из наиболее перспективных для интерпретации постановок является следующая: необходимо обратить полные магнитно-градиентные данные в трехмерной области с целью восстановления объемной намагниченности [13]. Подобная обратная задачи магнитостатики описывается операторным уравнением первого рода, которое в общем случае представляет собой некорректно поставленную задачу. Она может как иметь неединственное решение, так и не иметь классического решения вообще. При этом она, как правило, неустойчива по отношению к ошибкам измерения входных данных (компонент тензора ГКМИ). Эти трудности преодолеваются с помощью применения специальных методов решения таких некорректно поставленных задач – регуляризирующих алгоритмов (РА). Важнейшим классом регуляризирующих алгоритмов является семейство вариационных РА. Наиболее известен из них метод регуляризации Тихонова. Часто используется также обобщенный метод невязки.

Основная особенность постановки задачи, рассмотренной в [13], заключается в том, что рассматриваемая обратная задача является физически переопределенной. Классическая постановка задачи, заключающаяся в восстановлении магнитного момента тела по результатам измерения индуцированного им поля на некотором удалении от этого тела, является определенной, так как предполагает восстановление одной векторной функции по результатам измерений также одной векторной функции. Либо, принимая во внимание факт, что каждая компонента векторной функции является скалярной функцией, требуется восстановить три скалярные функции по результатам измерения трех скалярных функций (данная постановка приводит к системе из трех уравнений с тремя неизвестными функциями). Используя же постановку, которая связывает плотность магнитного момента тела с полным тензором ГКМИ, становится доступным дополнительное уравнение, которое связывает неизвестную векторную функцию с полным тензором градиентов компонент магнитной индукции, что, с учетом физических ограничений, дает нам дополнительные пять уравнений. В результате мы получаем физически переопределенную задачу, заключающуюся в определении 3 неизвестных скалярных функций по результатам экспериментальных измерений других 8 скалярных функций (подробности данной постановки рассмотрены в [13]).

Однако, если имеется априорная информация об индуцированном Землей магнитном поле в области расположения исследуемого в задаче тела, обратную задачу поиска векторной функции магнитного момента можно заменить на обратную задачу поиска скалярной функции магнитной восприимчивости. В этом случае физическая переопределенность модели будет более высокой, чем в модели, рассмотренной ранее [13]: будет требоваться найти 1 скалярную функцию по результатам экспериментальных измерений 8 других скалярных функций.

1. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

Уравнение, определяющее магнитное поле ${{{\mathbf{B}}}_{{{\text{field}}}}}$, индуцированное объектом с распределением намагниченности ${\mathbf{M}}(x,y,z)$ и локализацией в области $V$, имеет вид [14], [15]:

(1)
${{{\mathbf{B}}}_{{{\text{field}}}}}({{{\mathbf{r}}}_{s}}) = {{\nabla }_{s}}\iiint\limits_V {\left( {{{\nabla }_{s}}\frac{1}{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}},{\mathbf{M}}({\mathbf{r}})} \right)}\,d{v}.$
Здесь $\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right| = \sqrt {{{{(x - {{x}_{s}})}}^{2}} + {{{(y - {{y}_{s}})}}^{2}} + {{{(z - {{z}_{s}})}}^{2}}} $ – расстояние между точкой ${{{\mathbf{r}}}_{s}} = ({{x}_{s}},{{y}_{s}},{{z}_{s}})$, в которой расположен измеряющий поле ${{{\mathbf{B}}}_{{{\text{field}}}}}$ сенсор $s$, и точкой ${\mathbf{r}} = (x,y,z)$ области $V$, в которой расположен магнитный источник с суммарным магнитным моментом на единицу объема вещества ${\mathbf{M}}({\mathbf{r}})$; ${{\nabla }_{s}} \equiv \left( {\tfrac{\partial }{{\partial {{x}_{s}}}},\tfrac{\partial }{{\partial {{y}_{s}}}},\tfrac{\partial }{{\partial {{z}_{s}}}}} \right)$ – оператор вычисления градиента по переменным с индексом $s$.

Учтем отношение, связывающее намагниченность ${\mathbf{M}}({\mathbf{r}})$ области $V$ c магнитной восприимчивостью $\chi ({\mathbf{r}})$ [15], [17]:

(2)
${\mathbf{M}}({\mathbf{r}}) = \chi ({\mathbf{r}}){{{\mathbf{B}}}^{0}}({\mathbf{r}}),$
где ${{{\mathbf{B}}}^{0}}({\mathbf{r}})$ – индуцированное Землей магнитное поле.

Исходя из предположения, что в области $V$ индуцированное Землей магнитное поле ${{{\mathbf{B}}}^{0}}({\mathbf{r}})$ является постоянным, перепишем выражение (2) в виде

(3)
${\mathbf{M}}({\mathbf{r}}) = \chi ({\mathbf{r}})\left| {{{{\mathbf{B}}}^{0}}} \right|{\mathbf{l}},$
где ${\mathbf{l}} \equiv ({{l}_{x}},{{l}_{y}},{{l}_{z}})$ – единичный вектор, сонаправленный с постоянным индуцированным магнитным полем ${{{\mathbf{B}}}^{0}}({\mathbf{r}})$, $\left| {{{{\mathbf{B}}}^{0}}} \right|$ – величина этого поля.

Подставив (3) в (1), получим

(4)
${{{\mathbf{B}}}_{{{\text{field}}}}}({{{\mathbf{r}}}_{s}}) = {{\nabla }_{s}}\iiint\limits_V {\left( {{{\nabla }_{s}}\frac{1}{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}},\chi ({\mathbf{r}})\left| {{{{\mathbf{B}}}^{0}}} \right|{\mathbf{l}}} \right)d{v}}.$

Учтя, что $({{\nabla }_{s}}u,{{{\mathbf{l}}}_{s}}) \equiv \tfrac{{\partial u}}{{\partial {{{\mathbf{l}}}_{s}}}}$, перепишем (4) в виде

(5)
${{{\mathbf{B}}}_{{{\text{field}}}}}({{{\mathbf{r}}}_{s}}) = \left| {{{{\mathbf{B}}}^{0}}} \right|{{\nabla }_{s}}\iiint\limits_V {\frac{\partial }{{\partial {{{\mathbf{l}}}_{s}}}}\frac{1}{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}\chi ({\mathbf{r}})d{v}}.$

Поменяв местами операторы ${{\nabla }_{s}}$ и $\tfrac{\partial }{{\partial {{l}_{s}}}}$, перепишем (5) и получим

${{{\mathbf{B}}}_{{{\text{field}}}}}({{{\mathbf{r}}}_{s}}) = \left| {{{{\mathbf{B}}}^{0}}} \right|\iiint\limits_V {\frac{\partial }{{\partial {{{\mathbf{l}}}_{s}}}}{{\nabla }_{s}}\frac{1}{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}\chi ({\mathbf{r}})d{v}}.$

Учитывая, что

${{\nabla }_{s}}\frac{1}{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}} = \frac{{{\mathbf{r}} - {{{\mathbf{r}}}_{s}}}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{3}}}},$
получаем

${{{\mathbf{B}}}_{{{\text{field}}}}}({{{\mathbf{r}}}_{s}}) = \left| {{{{\mathbf{B}}}^{0}}} \right|\iiint\limits_V {\frac{\partial }{{\partial {{{\mathbf{l}}}_{s}}}}\frac{{{\mathbf{r}} - {{{\mathbf{r}}}_{{\mathbf{s}}}}}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{3}}}}\chi ({\mathbf{r}})d{v}}.$

Учитывая, что

$\frac{\partial }{{\partial {{{\mathbf{l}}}_{s}}}}\frac{{{\mathbf{r}} - {{{\mathbf{r}}}_{s}}}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{3}}}} = \frac{{3({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})({\mathbf{r}} - {{{\mathbf{r}}}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{\mathbf{l}}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{3}}}},$
получаем интегральное уравнение
${{{\mathbf{B}}}_{{{\text{field}}}}}({{{\mathbf{r}}}_{s}}) = \left| {{{{\mathbf{B}}}^{0}}} \right|\iiint\limits_V {{{{\mathbf{K}}}_{{{\text{field}}}}}({\mathbf{r}} - {{{\mathbf{r}}}_{s}};{\mathbf{l}})\chi ({\mathbf{r}})d{v}}$
с ядром ${{{\mathbf{K}}}_{{{\text{field}}}}}$, имеющим вид

${{{\mathbf{K}}}_{{{\text{field}}}}}({\mathbf{r}} - {{{\mathbf{r}}}_{s}};{\mathbf{l}}) = \frac{{3({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})({\mathbf{r}} - {{{\mathbf{r}}}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{\mathbf{l}}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{3}}}}.$

Преобразовав ${{{\mathbf{B}}}_{{{\text{field}}}}}({{{\mathbf{r}}}_{s}})$ к виду

$\begin{gathered} {{{\mathbf{B}}}_{{{\text{field}}}}}({{{\mathbf{r}}}_{s}}) = {{B}_{x}}({{{\mathbf{r}}}_{s}})i + {{B}_{y}}({{{\mathbf{r}}}_{s}})j + {{B}_{z}}({{{\mathbf{r}}}_{s}})k = \left( {\left| {{{{\mathbf{B}}}^{0}}} \right|\iiint\limits_V {\left( {\frac{{3({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})(x - {{x}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{{{l}_{x}}}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{3}}}}} \right)\chi ({\mathbf{r}})d{v}}} \right)i + \\ + \;\left( {\left| {{{{\mathbf{B}}}^{0}}} \right|\iiint\limits_V {\left( {\frac{{3({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})(y - {{y}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{{{l}_{y}}}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{3}}}}} \right){\kern 1pt} \chi ({\mathbf{r}})d{v}}} \right)j + \left( {\left| {{{{\mathbf{B}}}^{0}}} \right|\iiint\limits_V {\left( {\frac{{3({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})(z - {{z}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{{{l}_{z}}}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{3}}}}} \right){\kern 1pt} \chi ({\mathbf{r}})d{v}}} \right)k, \\ \end{gathered} $
и переобозначив переменные как $i = x,y,z$ и $p = ({{p}_{x}},{{p}_{y}},{{p}_{z}}) = ({{x}_{s}},{{y}_{s}},{{z}_{s}})$, мы можем записать выражения для компонент вектора ${{{\mathbf{B}}}_{{{\text{field}}}}}({{{\mathbf{r}}}_{s}})$ в следующем виде:

${{B}_{i}}({{{\mathbf{r}}}_{s}}) = \left| {{{{\mathbf{B}}}^{0}}} \right|\iiint\limits_V {\left( {\frac{{3({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})(i - {{p}_{i}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{{{l}_{i}}}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{3}}}}} \right)\chi ({\mathbf{r}})d{v}}.$

Беря производные от ${{B}_{i}}$ по пространственным переменным $i = x,y,z$ и $j = x,y,z \ne i$, мы можем получить следующие выражения для диагональных и недиагональных элементов тензора градиентов компонент магнитной индукции ${{{\mathbf{B}}}_{{{\text{tensor}}}}}$:

$\begin{gathered} {{B}_{{ii}}}({{{\mathbf{r}}}_{s}}) = \left| {{{{\mathbf{B}}}^{0}}} \right|\iiint\limits_V {\left( {\frac{{15({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}}){{{(i - {{p}_{i}})}}^{2}}}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{7}}}} - \frac{{6{{l}_{x}}(i - {{p}_{i}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{3({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}}} \right)\chi ({\mathbf{r}})d{v}}, \\ {{B}_{{ij}}}({{{\mathbf{r}}}_{s}}) = \left| {{{{\mathbf{B}}}^{0}}} \right|\iiint\limits_V {\left( {\frac{{15({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})(i - {{p}_{i}})(j - {{p}_{j}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{7}}}} - \frac{{3{{l}_{x}}(j - {{p}_{j}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{3{{l}_{y}}(i - {{p}_{i}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}}} \right)\chi ({\mathbf{r}})d{v}}. \\ \end{gathered} $

Обратим внимание на то, что мы определили полный тензор градиентов компонент магнитной индукции ${{{\mathbf{B}}}_{{{\text{tensor}}}}}$, который в отличие от индукции магнитного поля ${{{\mathbf{B}}}_{{{\text{field}}}}}$ (которая имеет 3 компоненты) имеет 9 компонент и может быть записан в следующей матричной форме:

${{{\mathbf{B}}}_{{{\text{tensor}}}}} = [{{B}_{{ij}}}]\left[ {\begin{array}{*{20}{c}} {\frac{{\partial {{B}_{x}}}}{{\partial x}}}&{\frac{{\partial {{B}_{x}}}}{{\partial y}}}&{\frac{{\partial {{B}_{x}}}}{{\partial z}}} \\ {\frac{{\partial {{B}_{y}}}}{{\partial x}}}&{\frac{{\partial {{B}_{y}}}}{{\partial y}}}&{\frac{{\partial {{B}_{y}}}}{{\partial z}}} \\ {\frac{{\partial {{B}_{z}}}}{{\partial x}}}&{\frac{{\partial {{B}_{z}}}}{{\partial y}}}&{\frac{{\partial {{B}_{z}}}}{{\partial z}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{B}_{{xx}}}}&{{{B}_{{xy}}}}&{{{B}_{{xz}}}} \\ {{{B}_{{yx}}}}&{{{B}_{{yy}}}}&{{{B}_{{yz}}}} \\ {{{B}_{{zx}}}}&{{{B}_{{zy}}}}&{{{B}_{{zz}}}} \end{array}} \right],$
где
$\tfrac{{\partial {{B}_{x}}}}{{\partial y}} = \tfrac{{\partial {{B}_{y}}}}{{\partial x}},\quad \tfrac{{\partial {{B}_{x}}}}{{\partial z}} = \tfrac{{\partial {{B}_{z}}}}{{\partial x}},\quad \tfrac{{\partial {{B}_{y}}}}{{\partial z}} = \tfrac{{\partial {{B}_{z}}}}{{\partial y}}\quad {\text{и}}\quad \tfrac{{\partial {{B}_{x}}}}{{\partial x}} + \tfrac{{\partial {{B}_{y}}}}{{\partial y}} + \tfrac{{\partial {{B}_{z}}}}{{\partial z}} = 0,$
что приводит к тому, что мы имеем 5 различных компонент матрицы тензора.

В результате мы приходим к системе из двух интегральных уравнений Фредгольма I рода:

(6)
$\begin{gathered} {{{\mathbf{B}}}_{{{\text{field}}}}}({{{\mathbf{r}}}_{s}}) = \iiint\limits_V {{{{\mathbf{K}}}_{{{\text{MI}}}}}\left( {{\mathbf{r}} - {{{\mathbf{r}}}_{s}};\left| {{{{\mathbf{B}}}^{0}}} \right|,{\mathbf{l}}} \right)\chi ({\mathbf{r}})d{v}}, \\ {{{\mathbf{B}}}_{{{\text{tensor}}}}}({{{\mathbf{r}}}_{s}}) = \iiint\limits_V {{{{\mathbf{K}}}_{{{\text{MGT}}}}}\left( {{\mathbf{r}} - {{{\mathbf{r}}}_{s}};\left| {{{{\mathbf{B}}}^{0}}} \right|,{\mathbf{l}}} \right)\chi ({\mathbf{r}})d{v}}, \\ \end{gathered} $
где ${{{\mathbf{B}}}_{{{\text{field}}}}} = {{[{{B}_{x}}{{B}_{y}}{{B}_{z}}]}^{{\text{T}}}}$ и ${{{\mathbf{B}}}_{{{\text{tensor}}}}} = {{[{{B}_{{xx}}}{{B}_{{xy}}}{{B}_{{xz}}}{{B}_{{yz}}}{{B}_{{zz}}}]}^{{\text{T}}}}$, а ядра интегральных уравнений ${{{\mathbf{K}}}_{{{\text{MI}}}}}$ и ${{{\mathbf{K}}}_{{{\text{MGT}}}}}$ имеют вид
${{{\mathbf{K}}}_{{{\text{MI}}}}}\left( {{\mathbf{r}} - {{{\mathbf{r}}}_{s}};\left| {{{{\mathbf{B}}}^{0}}} \right|,{\mathbf{l}}} \right) = \left| {{{{\mathbf{B}}}^{0}}} \right|\left[ {\begin{array}{*{20}{c}} {\frac{{3({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})(x - {{x}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{{{l}_{x}}}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{3}}}}} \\ {\frac{{3({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})(y - {{y}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{{{l}_{y}}}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{3}}}}} \\ {\frac{{3({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})(z - {{z}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{{{l}_{z}}}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{3}}}}} \end{array}} \right],$
и

${{{\mathbf{K}}}_{{{\text{MGT}}}}}\left( {{\mathbf{r}} - {{{\mathbf{r}}}_{s}};\left| {{{{\mathbf{B}}}^{0}}} \right|,{\mathbf{l}}} \right) = \left| {{{{\mathbf{B}}}^{0}}} \right|\left[ {\begin{array}{*{20}{c}} {\frac{{15({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}}){{{(x - {{x}_{s}})}}^{2}}}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{7}}}} - \frac{{6{{l}_{x}}(x - {{x}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{3({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}}} \\ {\frac{{15({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})(x - {{x}_{s}})(y - {{y}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{7}}}} - \frac{{3{{l}_{x}}(y - {{y}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{3{{l}_{y}}(x - {{x}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}}} \\ {\frac{{15({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})(x - {{x}_{s}})(z - {{z}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{7}}}} - \frac{{3{{l}_{z}}(x - {{x}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{3{{l}_{x}}(z - {{z}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}}} \\ {\frac{{15({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})(y - {{y}_{s}})(z - {{z}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{7}}}} - \frac{{3{{l}_{z}}(y - {{y}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{3{{l}_{y}}(z - {{z}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}}} \\ {\frac{{15({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}}){{{(z - {{z}_{s}})}}^{2}}}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{7}}}} - \frac{{6{{l}_{z}}(z - {{z}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}} - \frac{{3({\mathbf{l}},{\mathbf{r}} - {{{\mathbf{r}}}_{s}})}}{{{{{\left| {{\mathbf{r}} - {{{\mathbf{r}}}_{s}}} \right|}}^{5}}}}} \end{array}} \right].$

2. МЕТОД РЕШЕНИЯ

Если мы предположим, что $V \subset P = \{ (x,y,z){\kern 1pt} :\;{{L}_{x}}\;\leqslant \;x\;\leqslant \;{{R}_{x}},\;{{L}_{y}}\;\leqslant \;y\;\leqslant \;{{R}_{y}},\;{{L}_{z}}\;\leqslant \;z\;\leqslant \;{{R}_{z}}\} $ и система сенсорных плоскостей ограничена прямоугольным параллелепипедом $Q = \{ ({{x}_{s}},{{y}_{s}},{{z}_{s}}) \equiv $ $ \equiv (s,t,r){\kern 1pt} :\;{{L}_{s}}\;\leqslant \;s\;\leqslant \;{{R}_{s}},\;{{L}_{t}}\;\leqslant \;t\;\leqslant \;{{R}_{t}},\;{{L}_{r}}\;\leqslant \;r\;\leqslant \;{{R}_{r}}\} $, мы можем переписать систему (6) в следующем операторном виде:

(7)
${\mathbf{A}}\chi \equiv \int\limits_{{{L}_{x}}}^{{{R}_{x}}} \,\int\limits_{{{L}_{y}}}^{{{R}_{y}}} \,\int\limits_{{{L}_{z}}}^{{{R}_{z}}} \,{\mathbf{K}}(s,t,r,x,y,z)\chi (x,y,z)dxdydz = {\mathbf{B}}(s,t,r),$
где ${\mathbf{B}} = {{[{{B}_{x}}{{B}_{y}}{{B}_{z}}{{B}_{{xx}}}{{B}_{{xy}}}{{B}_{{xz}}}{{B}_{{yz}}}{{B}_{{zz}}}]}^{{\text{T}}}}$ и ${\mathbf{K}} = {{[{{{\mathbf{K}}}_{{{\text{MI}}}}}{{{\mathbf{K}}}_{{{\text{MGT}}}}}]}^{{\text{T}}}}$. В случае только магнитных данных: ${\mathbf{B}} = {{[{{B}_{x}}{{B}_{y}}{{B}_{z}}]}^{{\text{T}}}}$ и ${\mathbf{K}} = {{{\mathbf{K}}}_{{{\text{MI}}}}}$. В случае только градиентных данных: ${\mathbf{B}} = {{[{{B}_{{xx}}}{{B}_{{xy}}}{{B}_{{xz}}}{{B}_{{yz}}}{{B}_{{zz}}}]}^{{\text{T}}}}$ и ${\mathbf{K}} = {{{\mathbf{K}}}_{{{\text{MGT}}}}}$. В дальнейшем соответствующие модели будем называть как MI + MGT, MI и MGT-модели.

Далее будем предполагать, что $\chi \in W_{2}^{2}(P)$, $B \in {{L}_{2}}(Q)$, а оператор ${\mathbf{A}}$ с ядром ${\mathbf{K}}$ непрерывен и однозначен. Норма правой части уравнения (7) вводится следующим образом:

${{\left\| {\mathbf{B}} \right\|}_{{{{L}_{2}}}}} = \sqrt {\left\| {{{B}_{x}}} \right\|_{{{{L}_{2}}}}^{2} + \left\| {{{B}_{y}}} \right\|_{{{{L}_{2}}}}^{2} + \left\| {{{B}_{z}}} \right\|_{{{{L}_{2}}}}^{2} + \left\| {{{B}_{{xx}}}} \right\|_{{{{L}_{2}}}}^{2} + \left\| {{{B}_{{xy}}}} \right\|_{{{{L}_{2}}}}^{2} + \left\| {{{B}_{{xz}}}} \right\|_{{{{L}_{2}}}}^{2} + \left\| {{{B}_{{yz}}}} \right\|_{{{{L}_{2}}}}^{2} + \left\| {{{B}_{{zz}}}} \right\|_{{{{L}_{2}}}}^{2}} .$

Пусть вместо точно известных ${\mathbf{\bar {B}}}$ и оператора ${\mathbf{A}}$ известны их приближенные значения ${{{\mathbf{B}}}_{\delta }}$ и ${{{\mathbf{A}}}_{h}}$ такие, что ${{\left\| {{{{\mathbf{B}}}_{\delta }} - {\mathbf{\bar {B}}}} \right\|}_{{{{L}_{2}}}}}\;\leqslant \;\delta $, ${{\left\| {{\mathbf{A}} - {{{\mathbf{A}}}_{h}}} \right\|}_{{W_{2}^{2} \to {{L}_{2}}}}}\;\leqslant \;h$. При выписанных условиях задача является некорректной, для ее решения необходимо построить регуляризирующий алгоритм. Воспользуемся алгоритмом, основанным на минимизации функционала Тихонова

(8)
${{F}^{\alpha }}[\chi ] = \left\| {{{{\mathbf{A}}}_{h}}\chi - {{{\mathbf{B}}}_{\delta }}} \right\|_{{{{L}_{2}}}}^{2} + \alpha \left\| \chi \right\|_{{W_{2}^{2}}}^{2}.$
Для любого $\alpha > 0$ существует единственная экстремаль функционала Тихонова $\chi _{\eta }^{\alpha }$, $\eta = \{ \delta ,h\} $, реализующая минимум ${{F}^{\alpha }}[\chi ]$. Для выбора параметра регуляризации можно использовать алгоритм конечномерного обобщенного принципа невязки [18]. При выборе параметра $\alpha = \alpha (\eta )$ по обобщенному принципу невязки
$\rho (\alpha ) \equiv \left\| {{{{\mathbf{A}}}_{h}}\chi _{\eta }^{\alpha } - {{{\mathbf{B}}}_{\delta }}} \right\|_{{{{L}_{2}}}}^{2} - \mathop {\left( {\delta + h{{{\left\| {\chi _{\eta }^{\alpha }} \right\|}}_{{W_{2}^{2}}}}} \right)}\nolimits^2 = 0,$
$\chi _{\eta }^{\alpha }$ стремится при $\eta \to 0$ к точному решению задачи в норме $W_{2}^{2}$, а следовательно, и равномерно на $P$.

В качестве метода минимизации функционала Тихонова применяется метод сопряженных градиентов.

3. ЧИСЛЕННЫЕ АСПЕКТЫ АЛГОРИТМА

После введения конечно-разностной аппроксимации и сеток в областях $P$ и $Q$ приближенное решение $X$ (вектор, состоящий из сеточных значений функции $\chi $), которое реализует минимум функционала (8), может быть найдено как решение системы

(9)
$(A_{h}^{{\text{T}}}{{A}_{h}} + \alpha {{R}^{{\text{T}}}}R)X = A_{h}^{{\text{T}}}{{B}_{\delta }},$
где R – конечно-разностная аппроксимация оператора R: ${{\left\| \chi \right\|}_{{W_{2}^{2}}}} = {{\left\| {{\mathbf{R}}\chi } \right\|}_{{{{L}_{2}}}}}$, размерность матрицы A: $({{N}_{A}} \times N)$, матрицы R: $({{N}_{R}} \times N)$, вектора $X$: $(N \times 1)$.

Для решения системы (9) мы будем использовать реализацию метода сопряженных градиентов, предложенную в работе [20].

Пусть ${{X}^{{(s)}}}$ – минимизирующая последовательность, ${{p}^{{(s)}}}$, ${{q}^{{(s)}}}$ – вспомогательные векторы, ${{p}^{{(0)}}} = 0$, ${{X}^{{(1)}}}$ – любая допустимая точка. Тогда формулы метода сопряженных градиентов для поиска решения ${{X}^{{(N)}}}$ системы (9) могут быть записаны в виде

${{r}^{{(s)}}} = \left\{ \begin{gathered} A_{h}^{{\text{T}}}({{A}_{h}}{{X}^{{(s)}}} - {{B}_{\delta }}) + \alpha {{R}^{{\text{T}}}}(R{{M}^{{(s)}}}),\quad {\text{если}}\quad s = 1, \hfill \\ {{r}^{{(s - 1)}}} - {{q}^{{(s - 1)}}}{\text{/}}({{p}^{{(s - 1)}}},{{q}^{{(s - 1)}}}),\quad {\text{если}}\quad s\; \geqslant \;2, \hfill \\ \end{gathered} \right.$
${{p}^{{(s)}}} = {{p}^{{(s - 1)}}} + \frac{{{{r}^{{(s)}}}}}{{({{r}^{{(s)}}},{{r}^{{(s)}}})}},\quad {{q}^{{(s)}}} = A_{h}^{{\text{T}}}({{A}_{h}}{{p}^{{(s)}}}) + \alpha {{R}^{{\text{T}}}}(R{{p}^{{(s)}}}),\quad {{X}^{{(s + 1)}}} = {{X}^{{(s)}}} - \frac{{{{p}^{{(s)}}}}}{{({{p}^{{(s)}}},{{q}^{{(s)}}})}}.$

Необходимо отметить, что при проведении численных экспериментов мы определяем $\alpha = 0$, ${{X}^{{(1)}}} = 0$ и используем номер итерации $s$ в качестве параметра регуляризации. В этом случае критерий прекращения итерационного процесса согласуется с погрешностью задания входных данных посредством условия [21]

${{\left\| {{{A}_{h}}{{X}^{{(s + 1)}}} - {{B}_{\delta }}} \right\|}^{2}} \leqslant \mathop {\left( {\delta + h{{{\left\| {\chi _{\eta }^{\alpha }} \right\|}}_{{W_{2}^{2}}}}} \right)}\nolimits^2 .$

При численном решении системы (9) возможно эффективное применение многопроцессорных систем, особенности работы с которыми указаны в [4], [19].

4. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ ЭКСПЕРИМЕНТОВ

С целью сравнения различных моделей (MI, MI + MGT и MGT-модели) была проведена симуляция экспериментальных данных с относительной ошибкой $\delta \sim 4\% $ (как для MI-данных, так и для MGT-данных). Тестовые данные моделировались в области $P = \{ (x,y,z): - 5000\;\leqslant \;x\;\leqslant \;5000,$ $ - 5000\;\leqslant \;y\;\leqslant \;5000,$ $ - 105\;\leqslant \;z\;\leqslant \; - 95\} $ (что соответствует восстановлению распределения магнитной восприимчивости на глубине $z = - 100$) с сеткой (Nx, Ny, Nz) = (80, 80, 1) и области $Q = \{ ({{x}_{s}},{{y}_{s}},{{z}_{s}}) \equiv (s,t,r): - 4000\;\leqslant \;s\;\leqslant \;4000,$ $ - 4000\;\leqslant \;t\;\leqslant \;4000,$ ${{L}_{r}}\;\leqslant \;r\;\leqslant \;2000\} $ с сеткой (Ns, Nt, Nr) = (350, 20, 1). Модельная функция $\chi $ представлена на фиг. 1a. Результаты вычислений с использованием MI, MI + MGT и MGT-моделей представлены на фиг. 1б–г. Фиг. 1в демонстрирует, что MGT-модель дает результат с наилучшей детализацией решения.

Для тестовых вычислений использовалось оборудование Центра коллективного пользования сверхвысокопроизводительными вычислительными ресурсами МГУ им. М.В. Ломоносова [22], [23].

Фиг. 1.

Результаты вычислений: (a) модельное решение – функция $\chi (x,y)$, (б) восстановленное решение при использовании MI-модели, (в) восстановленное решение при использовании MGT-модели, (г) восстановленное решение при использовании MI + MGT-модели.

ЗАКЛЮЧЕНИЕ

Использование полного тензора градиентов компонент магнитной инукции (ГКМИ) для восстановления функции, определяющей магнитную восприимчивость, – важная задача в современной геофизике. Результаты численных расчетов продемонстрировали, что модель, использующая для решения поставленной обратной задачи полный тензор ГКМИ (MGT-модель), дает наилучший результат при восстановлении деталей модельного решения. Использование же полных магнито-градиентных данных (MI + MGT-модель) не дает никаких преимуществ в качестве восстановления решения по сравнению с использованием только магнитных данных (MI-модель).

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

  1. Lelievre P.G., Oldenburg D.W. Magnetic forward modelling and inversion for high susceptibility // Geophysical Journal International. 2006. V. 166. P. 76–90.

  2. Li Y.G., Oldenburg D.W. 3-D inversion of magnetic data // Geophysics. 1996. V. 61. P. 394–408.

  3. Pignatelli A., Nicolosi I., Chiappini M. An alternative 3D inversion method for magnetic anomalies with depth resolution // Annals of Geophysics. 2006. V. 49. P. 1021–1027.

  4. Lukyanenko D.V., Yagola A.G., Evdokimova N.A. Application of inversion methods in solving ill-posed problems for magnetic parameter identification of steel hull vessel // Journal of Inverse and Ill-Posed Problems. 2011. V. 18. № 9. P. 1013–1029.

  5. Christensen A., Rajagopalan S. The magnetic vector and gradient tensor in mineral and oil exploration // Preview. 2000. V. 84. P. 77.

  6. Heath P., Heinson G., Greenhalgh S. Some comments on potential field tensor data // Exploration Geophysics. 2003. V. 34. № 2. P. 57–62.

  7. Schiffler M., Queitsch M., Stolz R., Chwala A., Krech W., Meyer H.-G., Kukowski N. Calibration of SQUID vector magnetometers in full tensor gradiometry systems // Geophysical Journal International. 2014. V. 198. № 2. P. 954–964.

  8. Schmidt P.W., Clark D.A. Advantages of measuring the magnetic gradient tensor // Preview. 2000. V. 85. P. 26–30.

  9. Schmidt P.W., Clark D.A., Leslie K.E., Bick M., Tilbrook D.L. GETMAG-a SQUID magnetic tensor gradiometer for mineral and oil exploration // Exploration Geophysics. 2004. V. 35. P. 297–305.

  10. Zhdanov M.S., Cai H.Z., Wilson G.A. 3D inversion of SQUID magnetic tensor data // Geology and Geosciences. 2012. V. 1. P. 1–5.

  11. Wang Y., Rong L., Qiu L., Lukyanenko D.V., Yagola A.G. Magnetic susceptibility inversion method with full tensor gradient data using low-temperature SQUIDs // Petroleum Science. 2019. V. 16. № 4. P. 794–807.

  12. Ji S., Wang Y., Zou A. Regularizing inversion of susceptibility with projection onto convex set using full tensor magnetic gradient data // Inverse Problems in Science & Engineering. 2016. V. 4. P. 323–326.

  13. Wang Y., Lukyanenko D., Yagola A. Magnetic parameters inversion method with full tensor gradient data // Inverse Problems & Imaging. 2019. V. 13. № 4. P. 745–754.

  14. Zhdanov M.S. Integral transforms in geophysics // Springer Science & Business Media. 2012.

  15. Ji S., Wang Y., Zou A. Regularizing inversion of susceptibility with projection onto convex set using full tensor magnetic gradient data // Inverse Problems in Science & Engineering. 2017. V. 25. № 2. P. 202–217.

  16. Zhdanov M.S. Geophysical inverse theory and regularization problems // Elsevier. 2002. V. 36.

  17. Wang Y.F., Stepanova I.E., Titarenko V.N., Yagola A.G. Inverse problems in geophysics and solution methods // Higher Education Press, Beijing. 2011.

  18. Tikhonov A.N., Goncharsky A.V., Stepanov V.V., Yagola A.G. Numerical methods for the solution of ill-posed problems. Dordrecht: Kluwer Academic Publishers. 1995.

  19. Lukyanenko D.V., Yagola A.G. Some methods for solving of 3d inverse problem of magnetometry // Eurasian Journal of Mathematical and Computer Applications. 2016. V. 4. № 3. P. 4–14.

  20. Kalitkin N.N.,Kuz’mina L.V Improved form of the conjugate gradient method // Mathematical Models and Computer Simulations. 2012. V. 4. № 1. P. 68–81.

  21. Alifanov O.M., Artuhin E.A., Rumyantsev S.V. Extreme methods for the solution of ill-posed problems // Moscow: Nauka, 1988.

  22. Sadovnichy V., Tikhonravov A., Voevodin Vl., Opanasenko V. “Lomonosov”: Supercomputing at Moscow State University // In Contemporary High Performance Computing: From Petascale toward Exascale (Chapman & Hall/CRC Computational Science) Boca Raton, USA, CRC Press. 2013. P. 283–307.

  23. Voevodin Vl., Antonov A., Nikitenko D., Shvets P., Sobolev S., Sidorov I., Stefanov K., Voevodin Vad., Zhumatiy S. Supercomputer Lomonosov-2: Large Scale, Deep Monitoring and Fine Analytics for the User Community // Supercomputing Frontiers and Innovations. 2019. V. 6. № 2. P. 4–11.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал