Вариационно-комбинаторный подход к решению линейных и нелинейных обратных задач геофизики

Полный текст

ВВЕДЕНИЕ

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

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

Любой специалист в области теории интерпретации геофизических данных сталкивается с дилеммой: стремиться ли к условной сходимости решения конечно-разностных аналогов операторов Лапласа, Пуассона и т. п. к непрерывному аналогу при уменьшении шага сети наблюдений или строить некоторое множество оптимальных в определённом смысле решений изначально дискретных задач? Однозначного ответа на этот принципиально важный вопрос на наш взгляд нет. Всегда следует учитывать имеющуюся априорную информацию об элементах поля и его предполагаемых источниках. Даже спутниковые системы высокого разрешения, работающие практически в непрерывном режиме, не могут обеспечить исследователей такого рода дополнительными данными, которые гарантировали бы однозначную разрешимость трёхмерных обратных задач геофизики. В [5] В.Н. Страхов показал, что существует несчётное множество эквивалентных по внешнему гравитационному полю трёхмерных распределений масс, не отличающихся какой бы то ни было симметрией.

В связи с вышеизложенным, можно попытаться создать такой алгоритм действий при работе с большими и сверхбольшими геофизическими данными, который, с одной стороны, позволял бы максимально эффективно использовать всю полученную с помощью измерительной аппаратуры информацию об объектах и физических полях, а с другой стороны, включал бы в себя процедуры сортировки и анализа данных на основе научно обоснованного подхода. Поскольку мы имеем дело с дискретными наборами данных, то вариационные постановки желательно формулировать в сеточных пространствах и для сеточных же функций. Разнообразие вариантов расположения источников в узлах сетки, а также способов задания самих сеток может быть учтено путём определения некоторых новых функционалов качества сеточного решения – функционалов, значения которых на элементах векторного пространства (а таковыми являются сеточные функции и вектор-функции, заданные в узлах трёхмерных сеток) зависят от конфигурации узлов наблюдений в трёхмерном пространстве. Методы комбинаторной теории отличаются большим разнообразием: это и работа с графами, деревьями, лесами, и задание действия группы подстановок на элементах конечного множества, и различные представления объектов в специфических базисах и т.п. [6, 7]. Мы ограничимся в данной статье понятиями нумератора объектов некоторого дискретного множества и комбинаторной конфигурации.

1. ВАРИАЦИОННЫЕ ЗАДАЧИ В РАМКАХ ТЕОРИИ ДИСКРЕТНОГО ПОТЕНЦИАЛА С ПРИМЕНЕНИЕМ КОМБИНАТОРНОГО ПОДХОДА

Будем предполагать, что набор геофизических данных представляет собой дискретное множество $X_L=\left\{x_1,x_2\mathrm{,...,}{x}_L\right\}$, состоящее из L элементов. Элементами такого множества могут быть значения одного или нескольких физических полей в точках некоторой сети наблюдений. Ниже мы опишем более детально свойства элементов множества $X_L$.

Напомним читателю определение (b, v, r, k, λ)-конфигурации [6].

Система подмножеств X₁, Х₂,...,Х_b множества $X_v=\left\{x_1,x_2\mathrm{,...,}{x}_v\right\}$ называется (b, v, r, k, λ)-конфигурацией, если выполнены следующие условия:

1. $\left|X_i\right|=k,i=\mathrm{1,...,}b.$
2. Для любого $x_i$ существует ровно $r$ подмножеств $X_{{\mu }_1}...X_{{\mu }_r}$, которым принадлежит этот элемент.
3. Для любых $x_i,x_j,i\ne j,$ существует ровно $\lambda$ подмножеств $X_{v_1}...X_{v_{\lambda }}$, которым принадлежит пара элементов $x_i,x_j,i\ne j$.
4. Числа $v, k, \lambda$ удовлетворяют условиям $\lambda >0; k<v-1$.

С помощью так называемого метода “включения-исключения” [6] можно установить справедливость следующей теоремы.

Теорема 1 [6]. Пусть даны N элементов $a_1,a_2\mathrm{,...,}{a}_N$, которым приписаны некоторые веса $\omega \left(a_1\right), \omega \left(a_2\right),...,\omega \left(a_N\right)$, являющиеся элементами некоторого кольца К. Указанные элементы могут обладать или не обладать некоторыми свойствами $A_1, A_2, ..., A_N$. Обозначим суммарный вес элементов, обладающих ровно r свойствами, через M(r). Тогда верна следующая формула:

$M\left(r\right)=\sum _{k=r}^n{\left(-1\right)}^{k-r}\left(\begin{array}{c}k\\ r\end{array}\right)S_k, r=0,...,n.S_k=\sum _{1\le i_1\le 1_2...i_k\le n}M\left(A_{i_1}{A}_{i_2...}{A}_{i_k}\right)$, (1)

где $M\left(A_{i_1}{A}_{i_2}\mathrm{...}{A}_{i_k}\right)$ ‒ суммарный вес элементов, обладающих фиксированными свойствами $A_{i_1},A_{i_2}\mathrm{,...}{A}_{i_k}$ одновременно, $\left(\begin{array}{c}k\\ r\end{array}\right)$ ‒ число сочетаний из k по r. Формула (1) может быть полезна при представлении выборки измеренных значений некоторой физической величины как блок-схемы или (b, v, r, k, λ)-конфигурации. Ниже мы приведём примеры такого рода объектов изучения. В комбинаторике важное значение имеют функции выборок, называемые “нумераторами” [6].

Опишем нумератор для комбинаторной схемы [6] в некоммутативном симметричном базисе. Такого рода разбиение выборок данных встречается тогда, когда мы хотим разбить дискретное множество на некоторое число блоков.

Пусть X – дискретное множество, состоящее из m элементов. Обозначим через $T\left({\beta }_1, {\beta }_2,...,{\beta }_m\right)$ число разбиений множества X на непустые подмножества, имеющие вторичную спецификацию $||1^{{\beta }_1}{2}^{{\beta }_2}...m^{{\beta }_m}||, {\beta }_1+2{\beta }_2+m{\beta }_m=m$, т.е. содержащих ${\beta }_i$ подмножеств (блоков) размера i. Можно показать, что

$T\left({\beta }_1, {\beta }_2,...,{\beta }_m\right)=\frac{m!}{{\left(1!\right)}^{{\beta }_1}·{\left(2!\right)}^{{\beta }_2}...{\left(m!\right)}^{{\beta }_m}{\beta }_1!{\beta }_2!...{\beta }_N!}$; ${\beta }_1+2{\beta }_2+...+m{\beta }_m=m$. (2)

Если считать, что для построения адекватной реальности модели гравитационного или магнитного поля необходимо рассмотреть n зон ($n\le m$), в которых могут быть измерены значения поля, то при объёме выборки данных m (т.е. всего имеется m измерений, и эти измерения первоначально полагаются одинаково значимыми) число возможных размещений точек наблюдения будет равно [6] $T\left(m,n\right)$, а если просуммировать по всем $n\le m$, то получим следующее выражение:

${\omega }_m=\sum _{k=1}^{\infty }\frac{k^m}{2^{k-1}}, T\left(m,n\right)={∆}^n{0}^m, 0^m\equiv {∆}^n{x}^m\overline{){}_{x=0}}, ∆f\left(x\right)=f(x+1)-f\left(x\right)$

${\omega }_m=\sum _{n=0}^{m}T\left(m,n\right)$. (3)

В (3) через ${∆}^n$ обозначена n-я степень оператора приращения функции целочисленного аргумента при изменении аргумента на единицу.

Известна асимптотика для ${\omega }_m$:

${\omega }_m=\frac{1}{2}m!{\left({\log }_{2}e\right)}^{m+1}\left(1+o\left(1\right)\right),m\to \infty$. (4)

Теперь приведём формулу для нумератора элементов выборки при представлении в коммутативном несимметричном n-базисе. Мы выбрали такое представление потому, что оно наиболее адекватно, на наш взгляд, отражает свойства данных наблюдений, если изначально разбить всё окружающее планету пространство на зоны, удалённые от фиксированного центра на различные расстояния. В некоторых зонах точек наблюдения может не быть вовсе, в то время как в отдельных областях может быть сосредоточено достаточно большое число пунктов измерений физических величин. В [8] мы доказали ряд теорем единственности решения систем линейных алгебраических уравнений (СЛАУ), к которым редуцируются вариационные постановки в рамках метода локальных и региональных S-аппроксимаций. Однозначность определения решения СЛАУ является необходимым условием существования устойчивого к помехам во входных данных приближённого решения линейной обратной задачи геофизики (например, задачи определения плотностей простого и двойного слоёв, генерирующих такое же гравитационное или магнитное поле, как и истинные носители масс).

Пусть $\Lambda =\left({\Lambda }_1, {\Lambda }_2,...,{\Lambda }_n\right)$ ‒ упорядоченный набор последовательностей, ${\Lambda }_j\subseteq N_0=\left\{0,1,...\right\}$. Определим производящую функцию от n+1 переменной:

$\Phi (t,x_1,x_2,...x_n; \Lambda )=\prod _{j=1}^n\sum _{{\alpha }_j\in {\Lambda }_j}{t}^{{\alpha }_j}{x}_j^{{\alpha }_j}$, (5)

где суммирование распространяется на все ${\alpha }_j\in {\Lambda }_j$.

Если представить нумератор (5) в виде

$\Phi \left(t,x_1,x_2,...,x_n;\Lambda \right)=\sum _{m=0}^{\infty }{t}^m \sum _{{\alpha }_1+{\alpha }_2...+{\alpha }_m=m {\alpha }_j\in {\Lambda }_j} x_1^{{\alpha }_1}{x}_2^{{\alpha }_2}...x_n^{{\alpha }_n}$, (6)

где суммирование ведётся по всем ${\alpha }_1+{\alpha }_2+...+{\alpha }_n=m, {\alpha }_j\in {\Lambda }_j, {\alpha }_j\ge 0, {\alpha }_j\in \mathrm{\mathbb{Z}}$ — множеству целых чисел, то можно заключить, что m элементов некоторой выборки разбиваются на группы (максимум n групп), в каждой из которых содержится ${\alpha }_j$ предметов. Нумератор перечисляет такие группы, а коэффициент при $t^m$ задаёт “вес” конкретной группы в общей массе элементов.

В частном случае, когда ${\Lambda }_1={\Lambda }_2=...={\Lambda }_n=N_0$, мы можем получить следующее упрощённое выражение для (6) [6]:

$\Phi \left(t,N_0\right)=\sum _{m=0}^{\infty }{C}_{nm}{t}^m={\left(1-t\right)}^{-n}, C_{nm}=\left(\begin{array}{c}n+m-1\\ m\end{array}\right)$. (7)

В (7) $C_{nm}$ — это количество сочетаний выборок объёма m с n различными элементами с неограниченным числом повторений или, что то же самое, число размещений m тождественных частиц по n ячейкам без ограничения на число частиц в ячейке. Можно, разумеется, рассматривать и другие конструкции, но это мы предполагаем сделать в будущем. Смысл определения нумераторов (2) или (5)–(7) заключается в том, что мы можем представить себе степень неоднозначности решения обратной задачи, сформулированной в рамках некоторой аналитической модели физического поля. Мы не раз подчёркивали важность применения структурно-параметрического подхода при решении геофизических задач интерпретационного характера (см., например, [9]). Каждому источнику интенсивных аномалий гравитационного или магнитного полей соответствует при таком подходе свой блок в матрице СЛАУ и векторе правой части (последний также имеет блочный вид). Но самая большая проблема, возникающая при решении обратных задач геофизики, состоит в том, что заранее не известно, как подмножества исходной выборки связаны с геологическими структурами. Кроме того, разноточность и разнородность данных приводят к тому, что измеренные значения физических величин обладают, если можно так выразиться, неодинаковой ценностью для построения аналитической модели изучаемого поля (или полей, при совместной интерпретации информации, полученной с помощью различных геофизических методов исследования).

Следует подчеркнуть, что в (6) число групп n может быть как сравнимо с m, числом точек дискретного множества, так и быть намного меньше его.

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

Мы подчёркивали в наших предыдущих работах [8, 9] важность применения структурно-параметрического подхода при интерпретации данных различной природы. Можно сказать, что отдельные носители простого и двойного слоёв, с помощью которых аппроксимируются элементы потенциальных полей, соответствуют блокам, или группам, точек в некотором дискретном множестве Х. Каждому такому блоку соответствует свой парциальный вектор в общем векторе значений физической величины:

$X=\left({\mathbf{x}}_{{\alpha }_1}, {\mathbf{x}}_{{\mathrm{\alpha }}_2};...{\mathbf{x}}_{{\mathrm{\alpha }}_{\mathrm{n}}}\right), {\mathbf{x}}_{{\alpha }_k}\mathbf{=}\left({\mathit{x}}_{{\mathrm{i}}_{\mathrm{1}}}\mathit{,}\mathit{ }{\mathit{x}}_{{\mathit{i}}_{\mathit{2}}}\mathit{,}\mathit{.}\mathit{.}\mathit{.}\mathit{,}{\mathit{x}}_{{\mathit{i}}_{{\mathit{\alpha }}_{\mathit{k}}}}\right)\subset X$,

${\alpha }_1+{\alpha }_2+...+{\alpha }_n=m, {\alpha }_j\in {\Lambda }_j, {\alpha }_j\ge 0, {\alpha }_j\in \mathrm{\mathbb{Z}}$. (8)

Таким образом, каждый парциальный носитель простого или двойного слоёв в методе линейных интегральных представлений [4, 8] создаёт своё поле, которое затем вносит вклад в измеренные значения (мы считаем, что справедлив принцип суперпозиции физических полей: в случае сильных магнитных полей указанный принцип может нарушаться).

Перейдем теперь к описанию нового алгоритма интерпретации геофизических данных.

Как комбинаторные методы могут помочь при решении геофизических задач интерпретационного характера? Для начала, нужно определить свойства $A_1,A_2\mathrm{,...,}{A}_N$.

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

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

Алгоритм отбора данных

1. Определяем границы полигона.
2. Делим полигон на заданные подобласти в соответствии с характером задачи. Находим условный центр каждой подобласти и её диаметр.
3. Сортируем данные таким образом, чтобы в каждую выделенную на текущем этапе интерпретации подобласть попало некоторое (возможно, варьируемое) число значений физической величины.
4. В каждой подобласти находим 10 процентов элементов, имеющих минимальные и максимальные по модулю значения.
5. Формируем основные и контрольные массивы данных для каждой подобласти. В контрольные наборы должны входить как малые по модулю элементы, так и, наоборот, большие. Одновременно большие и малые значения в фиксированный контрольный набор не входят.
6. По основному набору данных для каждой подобласти строим математические модели физических полей. Качество аппроксимации проверяем на контрольных выборках для данной подобласти.
7. Формируем из векторов, принадлежащих к основным наборам, массив значений, по которому определяются аналитические аппроксимации изучаемой физической величины во всей области. В качестве контрольных точек для всей области берём 10 процентов минимальных по модулю элементов из всей выборки.
8. Сравниваем показатели качества решения для подобластей и всей области. “Выбрасываем” подобласти с плохими в определённом смысле показателями и снова строим аналитические аппроксимации для оставшейся части области. Сравниваем показатели качества решения для подобластей из этой части и всей части и т.д.
9. Процесс построения математической модели поля считается оконченным на данном этапе интерпретации в том случае, когда выбрасывание одной или нескольких подобластей не приводит к существенному (порядка 10‒15 процентов) улучшению качества решения интерпретационной задачи либо когда от всей области остаётся половина (последнее требование может быть изменено – это зависит от характера решаемой задачи).

Теперь попытаемся математически описать алгоритмический подход к решению обратных геофизических задач.

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

Пример № 1.

Если нам дано L значений некоторой физической величины, измеренной на некотором полигоне, то числа $(b, v, r, k, \lambda )$ будут иметь следующий смысл:

1. v – “отвечает” за расстояние от данной точки до выделенного фиксированного центра (например, центра подобласти);
2. r – точек находятся на одном и том же расстоянии при фиксированном v.
3. b – соответствует углу $\vartheta$ в сферической системе координат: у нас имеется столько различных значений этого угла;
4. k – точек лежат на конусе, характеризующемся фиксированным значением угла $\vartheta$; $vr=bk$;
5. Точки, соответствующие различным v, встречаются попарно ровно $\lambda$ на конусах.

Можно резюмировать: если точки наблюдения характеризуются сферическими координатами, то по выборке этих точек можно задать описанную выше (b, v, r, k, λ)-конфигурацию.

Пример №2.

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

Пусть b в данном случае – число некоторых опорных точек, или, центров исследования. В качестве свойств исследуемых величин, которыми они, эти величины, должны обладать примем: а) некоторый диапазон расстояний точки наблюдения от фиксированного центра; б) модуль вектора магнитной индукции (или некоторый диапазон модулей); в) модуль гравитационного поля (или также диапазон его изменения). Все точки наблюдения разобьются на блоки. Если различных расстояний от точки до какого-либо центра исследования имеется v, внутри s- го блока может находиться ровно k элементов, то всего точек, таким образом структурированных, будет $vr=bk=N_1\le N$. Каждое из v расстояний до центра исследования встречается ровно в r блоках.

Разбиение множества значений поля, описанное в данном примере, может оказаться полезным как в случае сферической системы координат, так и в декартовой. Важное отличие от примера № 1 состоит в том, что центры исследования не фиксированы, а определяются каждый раз заново, в зависимости от выборки или характера задачи. Если исходные данные связаны со сферической системой координат, а нам нужно изучить “поведение” физических полей вблизи нескольких (иногда – достаточно большого количества!) пунктов наблюдений, то целесообразно перейти в декартову (аффинную) систему координат.

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

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

При построении математических моделей полей наиболее “ценными” являются, безусловно, данные, полученные вблизи поверхности планеты. Поэтому такого рода данные целесообразно делить на группы элементов, которые отличаются угловыми координатами. Внутри каждой из групп можно задать нумератор вида (2), а затем найти наилучшее приближение к измеренному полю, строя аппроксимации элементов поля при различных конфигурациях узлов наблюдений, описываемых нумератором. Решение обратных задач следует находить, конечно, не при всех возможных разбиениях элементов на группы: число таких конфигураций при больших L огромно. Но какие-то наиболее важные для практических целей случаи необходимо рассмотреть. Речь идет об областях, подлежащих более детальному изучению: в них нужно взять больше элементов из выборки, чем в других зонах. Если же мы сформулируем вариационную постановку сразу для всего массива разнородных и разноточных данных (зачастую ещё и зависящих от времени), то качество приближённого решения в интересующей нас области может оказаться неудовлетворительным по причине неоднозначности решения обратной задачи.

1.1. ВАРИАЦИОННАЯ ПОСТАНОВКА ДЛЯ ОПРЕДЕЛЕНИЯ ГРАВИТАЦИОННОГО ПОЛЯ. ДЕКАРТОВА СИСТЕМА КООРДИНАТ

В [1–3] приводятся основные формулы теории дискретных гравитационного и магнитного потенциалов, а также примеры постановок вариационных задач по определению гравитационного и магнитного полей в сеточных пространствах.

Аналогами операторов Лапласа и Пуассона в декартовом сеточном пространстве $X_n^{\left(S\right)}$ (n измерений) являются следующие выражения:

$$\begin{gathered} \Lambda \left\{V_S(x^{\left(S\right)}\right\}=0; \\ \Lambda \left\{V_S(x^{\left(S\right)}\right\}=-\frac{G{C}_n{m}_s\left(x^{\left(S\right)}\right)}{H^{n-2}}, H={\left(\prod _{k=1}^n{h}_k\right)}^{\frac{1}{n}}, m_s\left(x^{\left(S\right)}\right)=H^n\rho \left(x^{\left(S\right)}\right). \end{gathered}$$(9)

В (9) $\Lambda$ ‒ выбранный конечноразностный аналог оператора Лапласа, а $m_s\left(x^{\left(S\right)}\right)=H^n\rho \left(x^{\left(S\right)}\right)$ — значения сеточных масс, порождающих сеточное гравитационное поле; $h_k$ — шаг по к-ому направлению в декартовой системе координат; $G$ — гравитационная постоянная; $C_n$ — константа, зависящая от размерности пространства, в трёхмерном пространстве: $C_3=4\mathrm{\pi }$

Оператор Лапласа в (9) будем считать заданным на шаблоне “крест” [2]. Тогда для трёхмерного сеточного пространства будем иметь:

$\begin{array}{l}{\Lambda }^{+}\left(V\right(x^{\left(S\right)}\left)\right)=-3\left(\sum _{p=1}^3{C}_p\right)V_{}^{\left(S\right)}\left(x^{\left(S\right)}\right)+\sum _{p=1}^3\left\{V_{}^{\left(S\right)}(x^{\left(S\right)}-h_p{e}_p)+V_{}^{\left(S\right)}(x^{\left(S\right)}+h_p{e}_p)\right\}=\mathrm{0,}\\ C_p=\frac{H^3}{h_1{h}_2{h}_3},H=\sqrt[3]{h_1{h}_2{h}_3},x^{\left(S\right)}=\left(x_1^{\left(s\right)},x_2^{\left(S\right)},x_3^{\left(S\right)}\right),x_p^{\left(S\right)}=k{h}_p,-K\le k\le +K;\\ h_p=const,p=\mathrm{1,2,3};\left|x_p^{\left(S\right)}\right|\le d_p=K{h}_p,p=\mathrm{1,2,3.}\end{array}$ (10)

В настоящей работе трёхмерное сеточное пространство (т.е. неограниченное множество точек $x^{\left(S\right)}={\left(x_1^{\left(S\right)}, x_2^{\left(S\right)}, ..., x_n^{\left(S\right)}\right)}^T$, $x_k^{\left(S\right)}=p{h}_k$; $-\infty <p<+\infty$, $p=\pm 0,\pm 1,\pm 2,...$; $h_k=const$, заменяется ансамблем, или семейством, расширяющихся компактов, представляющих собой прямоугольные параллелепипеды $\left|x_p^{\left(S\right)}\right|\le d_p=K_n{h}_p,p=\mathrm{1,2,3}$; здесь $d_p$, $K_n$ и $h_p$, $p=1,2,3$, ‒ некоторые положительные константы, причем $K_n\to \infty$ при $n\to \infty$.

Рассмотрим уравнения Лапласа и Пуассона в некотором прямоугольном параллелепипеде. Для того чтобы корректно поставить задачу по определению $V\left(x^{\left(S\right)}\right)$ в указанном параллелепипеде, нужно задать граничные условия. Будем считать, что на гранях параллелепипеда $K_{n_q}$ с номером $n_q$, принадлежащего описанному выше семейству расширяющихся компактов, значения гравитационного потенциала или его производных известны (континуальным аналогом дискретной постановки такого рода является задача Дирихле для уравнения Лапласа):

$V\left(x^{\left(S\right)}\right)=f\left(x^{\left(S\right)}\right)$. (11)

Решение задачи (10)‒(11) можно получить с помощью метода матричной прогонки. В указанном методе матрицы систем линейных алгебраических уравнений имеют трёхдиагональный или блочно-трёхдиагональный вид (если постановка рассматривается в трёхмерном сеточном пространстве).

Сформулируем теперь вариационную задачу для дискретного аналога оператора Лапласа в трёхмерном пространстве в декартовых координатах.

Дано. $D_S$ ‒ некоторая сеточная область в трёхмерном сеточном пространстве $X_3^{\left(S\right)}$; L точек наблюдения, в которых измерены значения гравитационного поля; ${\aleph }_{comb}\left(M_L\right)$ ‒ некоторое множество комбинаторных конфигураций, описываемых формулами (1) или (2), (5)–(7) (т.е. множество точек наблюдения M_L разбивается на подмножества, определяемые либо (b, v, r, k, λ)- конфигурацией, либо нумераторами (2), (5)–(7). Для ${\aleph }_{comb}\left(M_L\right)$ с помощью метода линейных интегральных представлений находятся аналитические аппроксимации гравитационного поля (и эквивалентные по внешнему полю распределения потенциалов простого и двойного слоёв на некотором семействе носителей) в соответствии с описанным выше алгоритмом отбора данных. Далее, определяются значения гравитационного поля на нескольких горизонтальных плоскостях ${\Pi }_k,k=\mathrm{1,...,}N$, расположенных над поверхностью планеты. Обозначим такое множество значений гравитационного поля через

$P_{grav}\left(\Pi \right)\equiv \left\{g\left(x^{\left(S\right)}\right),x^{\left(S\right)}\in {\Pi }_k,k=\mathrm{1,...,}N\right\},\Pi =\underset{k=1}{\overset{N}{\cup }}{\Pi }_k.$

Найти: минимум функционала $\Psi \left[V_S\left(x^{\left(S\right)}\right)\right]$, определяемого формулой

$\Psi \left[V_S\left(x^{\left(S\right)}\right)\right]={||{\Lambda }_{\Gamma }^{+}\left(V_S\left(x^{\left(S\right)}\right)\right)||}_{E\left(D_S\right)}^2+\alpha \Phi \left(P_{grav}^S\left(\Pi \right), V_S\left(x^{\left(S\right)}\right){|}_{\left(x^{\left(S\right)}\right)\in \Pi }\right), s\in {\aleph }_{comb}\left(M_L\right).$(12)

на множестве $V_S\left(x^{\left(S\right)}\right),x^{\left(S\right)}\in Int{D}_S$, где $\alpha \ge 0$ ‒ некоторый параметр (можно сказать, параметр регуляризации [9]); $E\left(D_S\right)$ ‒ эвклидова норма в сеточной области $D_S$; $\mathbf{\Phi }\left(t,z\right)$ ‒ неотрицательный функционал на векторах одинаковой размерности t и z. ${\Lambda }_{{}_{\Gamma }}^{+}$ ‒ это аппроксимация оператора Лапласа на шаблоне “крест” (см. (9)) вместе со значениями на гранях параллелепипеда $K_{n_q}$ с номером $n_q$, принадлежащего описанному выше семейству расширяющихся компактов, т.е. конечно-разностная аппроксимация оператора граничной задачи Дирихле для указанного параллелепипеда.

Таким образом, минимизация функционала (12) представляет собой поиск функции $V_S\left(x^{\left(S\right)}\right),x^{\left(S\right)}\in Int{D}_S$, такой чтобы она удовлетворяла уравнению Лапласа на шаблоне “крест” во внутренних узлах сеточной выбранной области, на границе выделенного параллелепипеда принимала бы заданные значения и, кроме того, совпадала бы со значениями проектора на каждую из горизонтальных плоскостей. При этом проектор на семейство горизонтальных плоскостей определяется для каждой конфигурации множества точек наблюдения, входящей в ${\aleph }_{comb}\left(M_L\right)$.

На практике, поиск минимума функционала (12) сводится к решению разреженной СЛАУ.

Важное замечание: условия на границе параллелепипеда $K_{n_q}$ можно также получить с помощью метода линейных интегральных представлений гравитационного поля (например, построив S-аппроксимации элемента поля [4] для заданной комбинаторной конфигурации множества точек наблюдения).

1.2. ВАРИАЦИОННЫЕ ПОСТАНОВКИ В СФЕРИЧЕСКОЙ СИСТЕМЕ КООРДИНАТ

Наряду с сеточными областями в трёхмерном евклидовом пространстве мы будем рассматривать их сферические аналоги:

$$\begin{gathered} {\mathbf{r}}^{\left(s\right)}\mathit{=}{\left(r^{\left(S\right)}, {\phi }^{\left(S\right)}, {\vartheta }^{\left(S\right)}\right)}^{\mathit{T}}\mathit{,}\mathit{ }{r}_{\mathit{k}}^{\mathit{\left(}\mathit{S}\mathit{\right)}}\mathit{=}\mathit{ }{r}_{\mathit{k}\mathit{-}\mathit{1}}^{\mathit{\left(}\mathit{S}\mathit{\right)}}\mathit{+}{h}_{\mathit{r}\mathit{,}\mathit{k}}\mathit{;}\mathit{ }k\mathit{=}\mathit{1}\mathit{,}\mathit{.}\mathit{.}\mathit{.}\mathit{;} \\ {\phi }_{\mathit{m}}^{\mathit{\left(}\mathit{S}\mathit{\right)}}\mathit{=}{\phi }_{\mathit{m}\mathit{-}\mathit{1}}^{\mathit{\left(}\mathit{S}\mathit{\right)}}\mathit{+}m{h}_{\mathit{\phi }}\mathit{,}\mathit{ }m\mathit{=}\mathit{1}\mathit{,}\mathit{.}\mathit{.}\mathit{.}\mathit{,}\mathit{ }{N}_{\mathit{2}}\mathit{,}\mathit{ }{h}_{\mathit{\phi }}\mathit{=}const\mathit{;} \\ {\vartheta }_{\mathit{l}}^{\mathit{\left(}\mathit{S}\mathit{\right)}}\mathit{=}{\vartheta }_{\mathit{l}\mathit{-}\mathit{1}}^{\mathit{\left(}\mathit{S}\mathit{\right)}}\mathit{+}{h}_{\mathit{\vartheta }\mathit{,}\mathit{l}}\mathit{,}\mathit{ }l\mathit{=}\mathit{1}\mathit{,}\mathit{.}\mathit{.}\mathit{.}\mathit{,}{N}_{\mathit{3}}\mathit{.} \end{gathered}$$ (13)

В (13) по радиусу и углу $\vartheta$ введены сетки с произвольными переменными шагами $h_{r,k}$ ,$h_{\vartheta ,i}$, и индексы k и l принимают значения из счётного и конечного множеств соответственно. Если сеточная область ограничена по радиусу и двум углам, будем полагать, что, $k=\mathrm{0,...,}{N}_1;l=\mathrm{0,...,}{N}_2;m=\mathrm{0,...,}{N}_3.$ Можно считать, что заданы три сеточные функции целочисленного аргумента $h_1\left(k\right)=h_{r,k}, h_2\left(m\right)=h_{\phi }, h_3\left(l\right)=h_{\vartheta ,l}$, значения которых в целых точках вещественной прямой равны шагам по соответствующим координатам.

Определим также сеточные функции $\rho \left(i\right)$ и $\tau \left(l\right)$ [10]:

$$\begin{gathered} \rho \left(k\right)=r_k, 1\le k\le N_1; \rho \left(0\right)=\left\{\begin{array}{l}\frac{1}{4}{h}_{r,1},r_0=0,\\ 0,r_0>0.\end{array}\right. \\ \tau \left(l\right)=\sin {\vartheta }_l, 1\le l\le N_3, \tau \left(0\right)=\left\{\begin{array}{l}\frac{1}{4}\sin h_{\vartheta ,1},{\vartheta }_0=0,\\ \sin {\vartheta }_0, \vartheta >0.\end{array}\right. \end{gathered}$$(14)

Запишем дискретный аналог оператора Лапласа в сферических координатах, действующего на сеточную функцию $V_{S,R}\left({}^{\left(S\right)}\right)$ (сеточная функция $V_{S,R}\left({}^{\left(S\right)}\right)$ может представлять собой сам потенциал либо какую-то его производную, умноженную на другую сеточную функцию, с тем чтобы в итоге произведение было гармоническим, например, $r^{\left(S\right)}{\partial }_{\hat{r}}\left(V_{S,R}\left({}^{\left(S\right)}\right)\right)$ ‒ гармоническая сеточная функция :

$$\begin{gathered} {∆}_h{V}_{S,R}\equiv \frac{1}{{\rho }^2} {\partial }_{\hat{r}}\left({\rho }^2{\partial }_{\overline{r}}{V}_{S,R}\right)+\frac{1}{{\rho }^2} {\partial }_{\hat{\vartheta }}\left(\tau \left(\vartheta \right) {\partial }_{\overline{\vartheta }} V_{S,R}\right)+\frac{1}{{\rho }^2\tau \left(\vartheta \right)} {\partial }_{\hat{\phi }}\left({\partial }_{\hat{\phi }}{V}_{S,R}\right); \\ h=\left(h_{r,k}, h_{\phi }, h_{\vartheta ,l}\right), k=1,...,N_1; l=1,...,N_3 \end{gathered}$$.(15)

В (15) по углу $\phi$ шаг считается постоянным и равным $h_{\phi }$. Сам угол $\phi$ принимает значения: ${\phi }_m=m{h}_{\phi }, m=1,...,N_2$. Разностные производные по координатам определяются следующим образом (приведём формулу для производной по радиусу):

${\partial }_{\overline{r}}\left(V_{S,R}\left(r^{\left(S\right)}\right)\right)=\frac{V_{S,R}(r_{k+1}^{\left(S\right)}, {\vartheta }_l^{\left(S\right)}, {\phi }_m^{\left(S\right)})-V_{S,R}(r_k^{\left(S\right)}, {\vartheta }_l^{\left(S\right)}, {\phi }_m^{\left(S\right)})}{h_{r,k}}$.

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

$$\begin{gathered} {\Psi }_R\left[V_{S,R}\left(r^{\left(S\right)}\right)\right]={||{\Lambda }_h\left(V_{S,R}\right(r^{\left(S\right)})||}_{E\left(D_{S,R}\right)}^2+\alpha \Phi \left(P_{grav}^S\left(\Xi \right),V_{S,R}(r^{\left(S\right)}{|}_{r^{\left(S\right)}\in \Xi }\right), \\ S\in {\aleph }_{comb}\left(M_L\right); \Xi ={\cup }_{k=1}^N{\Xi }_k \end{gathered}$$. (16)

В (16) $V_{S,R}\left(r^{\left(S\right)}\right)$ ‒ потенциал (или его производные!), записанный в сферических координатах, $\Phi \left(P_{{}_{grav}}^s\left(\Xi \right),{\left.V_{S,R}{(}^{\left(S\right)})\right|}_{{}^{\left(S\right)}\in \Xi }\right)$ ‒ билинейный функционал от своих аргументов; $P_{{}_{grav}}^s\left(\Xi \right)$ ‒ проектор на некоторое множество точек ; $\Xi =\underset{k=1}{\overset{N}{\cup }}{\Xi }_k$ — семейство концентрических сфер (можно задавать и другие геометрические структуры в сферических координатах, в точках которых вычисляются значения искомого элемента дискретного потенциала).

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

Предлагаемый в статье комбинированный подход к решению обратных линейных и нелинейных задач геофизики был апробирован на модельных и практических примерах. В модельном примере № 1 аномальное гравитационное поле создавалось набором из 10 прямых призм, нижние основания которых расположены на глубине 1 км под землёй, а высоты призм варьировались от 200 метров до 600 метров, а также 4 призм высотой 200 метров, нижние основания которых находились на глубине 3 км. Плотность породы в 6 верхних призмах полагалась равной 1.0 г/см³ и 0.5 г/см³ — в четырёх других; в нижних призмах плотность составила 1.0 г/см³. Участок съёмки имел размеры 100×100 км, рельеф местности выбирался спокойным, ровным (это означает, что перепад высот не превышал 250 метров). Пунктов наблюдения было 10 000 (т.е. вся выборка данных имела объем m = 10 000). Поперечники призм варьировались от 8 до 40 км. Поле призм аппроксимировалось суммой простого и двойного слоёв, распределённых на нескольких горизонтальных плоскостях (от 2 до 10), залегающих на глубинах от 100 метров до 500 м. Карта изолиний гравитационного поля призм приведена на рис. 1. Мы решили вариационную задачу (12) с помощью методики, описанной в [1] и нашли эквивалентные по внешнему гравитационному полю распределения масс на горизонтальных плоскостях для (b, v, r, k, λ)-конфигурации, в которой b = 10, k = 1 000, v = 20, r = 500, λ = 5. Поясним смысл параметров указанной конфигурации, которая соответствует примеру № 2 из первого раздела статьи.

 

Рис. 1. Гравитационное поле в модельном примере № 1. N = 10000.

 

b равно числу центров исследования, а именно центров верхних призм — они располагались на диагонали (см. рис. 1), с интервалом в 20 км; r —число блоков (центров исследования), в которых содержались точки, расположенные от центра на заданных расстояниях (от 5 до 20 км); v — число различных расстояний, на которых точки могут находиться в зоне, окружающей центр исследования; k — число точек в одном блоке; λ — число блоков, в которых могут находиться пары различных элементов. Мы получили ожидаемый результат: на показатель качества решения задачи по построению аналитических аппроксимаций гравитационного поля оказывали влияние только те точки призм, которые находились вблизи пунктов наблюдения. Показатель качества решения составил $∆={10}^{-10}$:

$∆=\frac{{||Az-f_{\delta }||}_E}{{||f_{\delta }||}_E}$,

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

Мы последовательно исключали из выборки блоки, соответствующие различным центрам исследования. Выяснилось, что показатель качества решения незначительно ухудшался (менее, чем в 10 раз) при отбрасывании блоков с меньшей по абсолютному значению плотностью породы и небольших по размеру. Если удалить три верхних блока (см. рис. 1), то показатель качества решения возрастёт в 10 000 раз. На рис. 2 изображено гравитационное поле, создаваемое двумя группами призм в модельном примере № 2. Было сгенерировано восемь призм в верхнем эшелоне и три — в нижнем. Плотности верхних призм чередовались: четыре призмы были заполнены веществом с плотностью 1.0 г/см³, в четырёх оставшихся плотность была равна 0.5 г/см³. В трёх нижних призмах вещество имело плотность 1.0 г/см³. Результаты расчётов в этом примере получились следующие: $∆={10}^{-12}$ для такой же комбинаторной конфигурации, что и в модельном примере № 1.

 

Рис. 2. Гравитационное поле в модельном примере № 2. N = 10000.

 

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

Данные миссии Messenger [11] за несколько дней полета в 2014 г. были представлены как некоторые комбинаторные конфигурации. Всего точек в наборе было 14 000. В файлах, содержащих “сырые данные”, указывались декартовы координаты точек наблюдения в километрах, при этом начало системы координат совпадает с центром масс Меркурия. Носители простого и двойного слоёв при аппроксимации лапласовых полей располагались в коре Меркурия, т.е. на расстоянии от 0.1 до 100 км от поверхности планеты; аппроксимировали мы компоненту $B_z$.

Мы описали в [12] алгоритм построения различных вариантов S-аппроксимаций компонент магнитного поля Меркурия. При работе с “сырыми” данными выяснилось, что те точки орбиты станции Messenger, которые находились на большом удалении от центра планеты (более 1000–1500 км), практически не влияли на качество аппроксимации: элементы матрицы СЛАУ, соответствующие таким пунктам наблюдений, имеют во много раз меньший модуль, чем остальные. Необходимы дальнейшие исследования в области изучения зависимости свойств орбиты спутника (эксцентриситета, возможных отклонений от эллиптичности и т.п.) и характеристических чисел матриц СЛАУ, к решению которых редуцируются задачи по построению аналитических моделей полей и нахождению эквивалентных по этим полям распределений источников.

Мы отметили в предыдущем разделе статьи, что при минимизации функционалов (12) и (16) можно применять метод матричной прогонки. Однако эквивалентные указанным вариационным постановкам СЛАУ достаточно хорошо решаются с помощью блочного метода регуляризации Холецкого (БМХР) и усовершенствованного блочного метода решения СЛАУ (УБМ) [13]. При этом мы полагали, что Меркурий представляет собой шар радиуса $R_0=2439.0$ км. Результаты аппроксимации представлены в табл. 1.

 

Таблица 1. Модифицированные S- аппроксимации z-компонентs магнитного поля Меркурия по данным Messenger в рамках комбинаторного подхода

№	$\frac{B}{N}$	$R_0$, км или $H_1, H_2$	Метод решения СЛАУ	${\mathrm{S}}_{\min }$, нТл	${\mathrm{S}}_{\max }$, нТл	${\mathrm{S}}_0$, нТл	$\frac{∆}{t}$
1	$$\begin{gathered} B_z \\ 14000 \end{gathered}$$	2400	БМХР	0.012	0.024	0.062	$$\begin{gathered} 2.3\times {10}^{-2} \\ 43:27 \end{gathered}$$
2	$$\begin{gathered} B_z \\ 5000, 5000, 2000 \end{gathered}$$	2400	БМХР	0.0027	0.0039	0.0033	$$\begin{gathered} 1.36\times {10}^{-6} \\ 26:31 \end{gathered}$$
3	$$\begin{gathered} B_z \\ 3000, 3000, 3000, 1000 \end{gathered}$$	2400	УБМ	0.001	0.0015	0.0012	$$\begin{gathered} 1.32\times {10}^{-8} \\ 27:13 \end{gathered}$$
4	$$\begin{gathered} B_z \\ 2000, 2000, 3000 \end{gathered}$$	2400	УБМ	0.001	0.0015	0.0011	$$\begin{gathered} 1.40\times {10}^{-8} \\ 32:18 \end{gathered}$$
5	$$\begin{gathered} B_z \\ 7000 \end{gathered}$$	2400	УБМ	0.01	0.015	0.012	$$\begin{gathered} 1.21\times {10}^{-4} \\ 20:37 \end{gathered}$$
6	$$\begin{gathered} B_z \\ 2000\times 7 \end{gathered}$$	2400	УБМ	0.001	0.0015	0.0014	$$\begin{gathered} 1.1\times {10}^{-7} \\ 33:42 \end{gathered}$$
Примечание. Обозначения, принятые в таблице 1. Здесь $\sigma =\frac{{||Ax-f_{\delta }||}_E}{\sqrt{N}}$ – среднеквадратическое отклонение, ${\sigma }_{max}=\sqrt{\frac{{\delta }_{max}^2}{N}}$, ${\sigma }_{min}=\sqrt{\frac{{\delta }_{min}^2}{N}}$, $∆=\frac{{||Ax-f_{\delta }||}_E}{{||f_{\delta }||}_E}$ – показатель качества решения, ${\sigma }_0$ - среднеквадратическое отклонение, полученное в результате решения СЛАУ, $t$ – время в часах, минутах и секундах.

 

Орбита космической миссии представляла собой достаточно вытянутый эллипс (самые удалённые точки траектории Messenger располагались на расстоянии порядка 20000 км). Поэтому для построения аппроксимаций мы отбирали только те точки, которые были от поверхности планеты не далее, чем в 800 км. Для выделения из “сырых” данных составляющих магнитной индукции, генерируемых токами в жидком ядре и коре – так называемого внутреннего магнитного поля Меркурия – можно воспользоваться приближением the thin shell approximation [14]. При таком подходе точки наблюдения должны находиться в пределах тонкой (по сравнению с некоторыми параметрами, характеризующими топологию планеты) оболочки, окружающей Меркурий. Полоидальное и тороидальное магнитные поля, создаваемые токами в плазме вокруг Меркурия, не вносят вклад в наблюдаемое поле в этом случае. Для того чтобы в наборе, по которому выполнялись аппроксимации, оставалось 11 000–14 000 элементов, мы синтезировали магнитное поле в промежуточных узлах сети наблюдений. По найденным распределениям эквивалентных источников мы находили пространственное распределение элементов магнитного поля,

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

Весь набор данных по магнитному полю Меркурия рассматривался нами как некоторое дискретное множество. Это дискретное множество было затем представлено шестью различными комбинаторными конфигурациями (они соответствуют строкам табл. 1, начиная с самой верхней): 1-я конфигурация — это вся выборка; вторая конфигурация — это разбиение множества точек на три блока (см. формулу (2)) указанных в таблице размеров; третья конфигурация — это разбиение выборки на 4 блока по формуле (2); четвёртая — это конфигурация, получающаяся при группировании точек наблюдения по зонам, удалённым от Меркурия на расстояния 10, 50 и 100 км; пятая конфигурация — это точки, расположенные на расстоянии 50 км от поверхности Меркурия; шестая конфигурация — это группы точек, располагающиеся в семи секторах, биссектрисы которых изображены чёрным цветом на рис. 1.

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

 

Рис. 3. Магнитное поле Меркурия по данным Messenger.

 

Рис. 4. Распределение эквивалентных по магнитному полю диполей на сфере R = 2400, построенное по всей выборке.

 

Рис. 5. Распределение эквивалентных по магнитному полю диполей на сфере R = 2400, построенное по аналитически построенному магнитному полю в точках орбиты, расположенных в первых трёх зонах (по степени удалённости от поверхности планеты): от 0 до 10 км; от 10 до 50 км и от 50 до 100 км.

 

ЗАКЛЮЧЕНИЕ

1. В статье рассмотрен новый подход к решению обратных линейных и нелинейных задач геофизики, ключевой особенностью которого является предварительное задание комбинаторных схем на основе исходной выборки и последующая формулировка вариационных задач в рамках теории дискретного потенциала и метода линейных интегральных представлений. Таким образом, построение математической модели физического поля подразумевает не просто обработку “сырых” геофизических данных и решение некоторой условно-вариационной задачи, а создание адекватного реальности алгоритма оценки и отбора данных с поэтапным анализом результатов решения “парциальных”, или вспомогательных, задач.
2. Дополнительные стабилизаторы, определяемые в статье, позволяют повысить показатель качества аналитической аппроксимации элемента гравитационного или магнитного полей. Проектирование элементов поля, восстановленных с помощью известной методики локальных и региональных модифицированных S-аппроксимаций, в точки некоторой системы подмножеств полупространства, свободного от источников поля, обеспечивает устойчивость приближённого решения обратной задачи по нахождению эквивалентных носителей.
3. В статье подчёркивается необходимость рассмотрения ансамбля конфигураций точек наблюдений и выявления априорных связей между значениями изучаемого объекта (поля) в различных областях как реального пространства $R^3$, так и его сеточного аналога.

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

ИСТОЧНИК ФИНАНСИРОВАНИЯ

Работа выполнена в рамках госзадания ИФЗ РАН.

Об авторах

И. Э. Степанова

Институт физики Земли им. О.Ю. Шмидта Российской Академии наук

Автор, ответственный за переписку.
Email: tet@ifz.ru
Россия, Москва

И. И. Колотов

Московский государственный университет им. М.В. Ломоносова

Email: tet@ifz.ru

физический факультет

Россия, Москва

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

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