Журнал вычислительной математики и математической физики, 2023, T. 63, № 8, стр. 1317-1331

ВВЕДЕНИЕ

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

1. ТЕОРЕМЫ ЕДИНСТВЕННОСТИ В СЛУЧАЕ ВОССТАНОВЛЕНИЯ МАГНИТНОЙ ИНДУКЦИИ

Приведем сначала выражения для магнитного поля (см. [2], [3]) в предположении, что среда является парамагнитной и измерения выполняются в некоторой декартовой системе координат. Если рассматриваются небольшие по площади полигоны (территории, на которых проводятся геофизические исследования), то можно считать поверхность Земли плоской, а саму планету – нижним полупространством. Постановки обратных задач при таких допущениях будем в дальнейшем называть локальными. Сферичность Земли обычно не учитывается, если линейный размер полигона не превышает 111 км, что соответствует приблизительно одному градусу по широте:

Здесь ${{\bar {r}}_{s}}$ – радиус-вектор точки наблюдения ${{P}_{s}}$ в декартовой системе координат, $\bar {r}$– радиус-вектор точки внутри носителя V, $\bar {B}({{\bar {r}}_{s}})$ – магнитная индукция в системе СИ, измеренная в точке, характеризующейся радиус-вектором ${{\bar {r}}_{s}}$, ${{\mu }_{0}}$ – магнитная постоянная в системе СИ. Через $K({{\bar {r}}_{s}},\bar {r})$ обозначена матричная функция координат точки наблюдения и точки, принадлежащей носителю, $M(\bar {r})$ – вектор-функция, характеризующая намагниченность среды.

Поле, создаваемое одним магнитным диполем, выражается формулой

Смысл обозначений в (1.2) такой же, как и в (1.1). Через $\bar {m} = {{\left( {{{m}_{1}},{{m}_{2}},{{m}_{3}}} \right)}^{{\text{т}}}}$ обозначен магнитный момент точечного диполя. Будем считать, что в точках наблюдения установлены условные “сенсоры”. Возникает вопрос об однозначной разрешимости следующей системы линейных алгебраических уравнений, возникающей, как уже отмечалось выше, в процессе дискретизации задачи (1.1):

В (1.3) через ${{\rho }_{{ki}}}$ обозначено расстояние между k-м сенсором и i-м диполем:

${{\bar {m}}_{i}} = {{\left( {{{m}_{{ix}}},{{m}_{{iy}}},{{m}_{{iz}}}} \right)}^{{\text{т}}}}$ – магнитный момент i-го диполя в декартовой системе координат. Можно утверждать, что система (1.3) имеет блочный вид где блок ${{K}_{{il}}}$ соответствует i-му сенсору и l-му диполю, ${{\bar {m}}_{i}}$ – вектор компонент i-го диполя, ${{\bar {B}}_{l}}$ – вектор-функция, измеренная в l-м сенсоре. Блок K_il имеет следующий вид:

Рассмотрим сначала случай двух сенсоров и двух магнитных диполей.

2. ТЕОРЕМЫ ЕДИНСТВЕННОСТИ В МЕТОДЕ S-АППРОКСИМАЦИЙ. ЛОКАЛЬНЫЙ ВАРИАНТ

Если известны компоненты магнитного или гравитационного полей (например, первая производная потенциала по z на некотором рельефе), то можно представить потенциал поля в виде суммы простого и двойного слоев, создаваемых совокупностью горизонтальных плоскостей, расположенных ниже заданного рельефа. Если систему координат выбрать так, чтобы дневная поверхность (поверхность Земли в локальном варианте метода S-аппроксимаций, см. [4], [5], [9–13]) задавалась уравнением z = 0, то потенциал представляется в виде (см. [8])

Тогда производная по z потенциала V, взятая с обратным знаком, будет иметь вид

Функции ${{\rho }_{{1,l}}},\;{{\rho }_{{2,l}}}$ неизвестны. Пусть компоненты поля заданы в конечном множестве точек М_i, M_i = $\;({{x}_{i}},{{y}_{i}},{{z}_{i}})$, i = 1, …, N. Обозначим подынтегральную функцию в первом слагаемом для l-го слоя в (2.2) в точке М_i через $Q_{{1,l}}^{{(i)}}$, а во втором слагаемом – через $Q_{{2,l}}^{{(i)}}$. Тогда получим

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

получим, что искомые функции должны иметь вид (см. [4], [5], [9])

Таким образом, приходим к следующей системе линейных уравнений:

элементы матрицы которой в нашем случае имеют вид

Элементы а_ij матрицы А могут быть вычислены явно с помощью интеграла Пуассона:

Поясним, как от формулы (2.8) перейти к выражению (2.9). На каждой из плоскостей${{H}_{r}}$ распределены источники гравитационного поля. Первое слагаемое в сумме (2.2) при фиксированном значении индекса l представляет собой двойной слой с плотностью ${{\rho }_{{1,l}}}(\xi )$, а второе слагаемое – производную этого двойного слоя по z, плотность распределения этой производной равна ${{\rho }_{{2,l}}}(\xi )$. Далее, в формуле (2.8) при каждом значении индекса l выражение

есть произведение ядра интеграла Пуассона, вычисленного в точке M_i с координатами $\;({{x}_{i}},{{y}_{i}},{{z}_{i}})$ и точке плоскости $\hat {\xi } = ({{\xi }_{1}},{{\xi }_{2}})$, на значение плотности ${{\rho }_{{1,l}}}(\xi )$, которая зависит от координат точки ${{M}_{j}} = ({{x}_{j}},{{y}_{j}},{{z}_{j}})$, следующее:

Таким образом, на плоскости ${{H}_{r}}$ значение искомой гармонической в верхнем полупространстве функции равно ${{\rho }_{{1,l}}}(\hat {\xi })$ в точке с координатами $\hat {\xi } = ({{\xi }_{1}},{{\xi }_{2}})$. Если теперь “подняться выше”, в точку $\;{{M}_{i}} = ({{x}_{i}},{{y}_{i}},{{z}_{i}})$, то к числителю дроби в выражении для $\rho _{{1,l}}^{{}}(\hat {\xi })$ нужно прибавить ${{z}_{j}} - {{H}_{l}}$, и на такое же значение изменится z-координата в знаменателе. Аналогичным образом трактуется выражение $Q_{{2,l}}^{{(i)}}(\hat {\xi })Q_{{2,l}}^{{(j)}}(\hat {\xi })$ в (2.8). Разница заключается только в том, что одна производная “перебрасывается” с ядра Пуассона на функцию, описывающую плотность масс в интеграле. Поэтому в (2.9) при каждом значении индекса, по которому выполняется суммирование, складываются первые и третьи производные потенциалов точечных источников.

Иначе говоря, в (2.9) при фиксированном значении индекса r стоит сумма первых и третьих производных по z от потенциала точечного источника, расположенного в точке ${{\tilde {M}}_{j}} = ({{x}_{j}},{{y}_{j}},2{{H}_{l}} - {{z}_{j}})$, вычисленных в точке $\;{{M}_{i}} = ({{x}_{i}},{{y}_{i}},{{z}_{i}})$.

По найденным из решения системы (2.7)–(2.9) множителям λ_i, I = 1, …, N, можно далее определить величины некоторых функционалов р_s, s = 1, …, S (см. [4], [9–13]). Эти функционалы могут представлять собой высшие производные потенциала, аналитическое продолжение гравитационного поля и т.д. Для приложений очень важно иметь возможность пересчета значений поля на некоторую регулярную сеть наблюдений.

Рассмотрим случай, когда носитель потенциалов простого и двойного слоев только один. Пусть D – плоскость, описываемая уравнением

По найденным из решения системы (2.7)–(2.9) плотностям ${{\rho }_{1}}({{\xi }_{1}},{{\xi }_{2}}) \in {{L}_{2}}(D)$ и ${{\rho }_{2}}({{\xi }_{1}},{{\xi }_{2}}) \in {{L}_{2}}(D)$ простого и двойного слоев, распределенных на этой плоскости, мы можем найти совокупность ограниченных линейных функционалов вида

где для всех s

Действительно, значения аномалий силы тяжести $g({{M}_{i}};\rho )$выражаются соотношением (2.10) при s = i:

а значения элементов аномального поля $u({{M}_{i}};\rho )$ – соотношением вида где $\alpha = ({{\alpha }_{1}},{{\alpha }_{2}},{{\alpha }_{3}})$ – мультииндекс, ${\text{|}}\alpha {\text{|}} = {{\alpha }_{1}} + {{\alpha }_{2}} + {{\alpha }_{3}} \geqslant 1,$и соотношением (2.10) при s = i:

Везде $\gamma $ – гравитационная постоянная ($\gamma = 6.67 \times {{10}^{{ - 11}}}$ Н ⋅ м²/кг²). Здесь точка ${{M}_{i}}$ имеет координаты ${{M}_{i}} = ({{x}_{i}},{{y}_{i}},{{z}_{i}}),$ ${{{\mathbf{r}}}_{{\mathbf{i}}}} = \left( {{{x}_{i}},{{y}_{i}},{{z}_{i}}} \right)$. Причем $z \geqslant - H$, т.е. мы вычисляем значения производных потенциала выше носителя, которым в нашем случае является плоскость D. Здесь $R({\mathbf{\xi }} - {{{\mathbf{r}}}_{i}})$ – длина вектора $({\mathbf{\xi }} - {{{\mathbf{r}}}_{i}})$.

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

Если какая-либо высшая производная гравитационного потенциала (см. (2.10)–(2.13)) аппроксимируется потенциалом простого слоя, то матрица системы (2.7)–(2.9) принимает вид

Попытаемся выяснить, по аналогии с разд. 1 настоящей статьи, в каких случаях матрица (2.16) является вырожденной. Для того чтобы определитель матрицы был равен нулю, достаточно, чтобы какие-либо две строки этой матрицы были пропорциональны друг другу, например, i-я и j-я. Рассмотрим более подробно элементы указанных строк и запишем условие пропорциональности:

Из (2.18) следует, что $r_{{ij}}^{2}$ всегда меньше или равно нулю:

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

Таким образом, мы можем сделать вывод о том, что верна

Теорема 5. Матрица системы (2.7)–(2.9) при представлении поля с помощью потенциала простого слоя имеет ранг не меньше двух; таким образом, для двух различных точек наблюдения система является невырожденной.

Перейдем теперь к случаю б) представления элементов гравитационного и магнитного полей в виде потенциала двойного слоя.

Матрица системы (2.9) приобретает следующий вид:

Запишем, как и для случая простого слоя, условие пропорциональности диагональных и внедиагональных элементов каких-нибудь двух строк:

(мы рассматриваем константу С как функцию координаты: $C = C({{z}_{j}})$).

Явную зависимость $r_{{ij}}^{2} = f(R,C({{z}_{j}}))$ найти трудно, поэтому попробуем определить, выполняется ли условие (2.20) на какой-либо сфере радиуса $R = \sqrt {{{{({{z}_{i}} + {{z}_{j}} + 2H)}}^{2}} + r_{{ij}}^{2}} $. Если некоторая точка трехмерного пространства лежит на такой сфере, то мы можем выразить ее координату z:

Будем считать, что ${{z}_{i}},\;{{z}_{j}}$ известны. Тогда ${{С}_{1}} \equiv {{z}_{i}} + {{z}_{j}} + 2H = {{C}_{1}}(C),$ $C = \frac{6}{{16{{{({{z}_{i}} + H)}}^{2}}{{{({{z}_{j}} + H)}}^{2}}}}$.

Исследуем, всегда ли существует такое $r_{{ij}}^{2}$ > 0 (при всех ли ${{C}_{1}}(C)$ и С), что справедливо равенство

Попробуем изобразить поиск решения графически (см. фиг. 4).

Если $С_{1}^{2} + r_{{ij}}^{2} = {{R}^{2}}$, то, согласно (2.20), на сфере радиуса R мы будем иметь:

Но $R \geqslant {{C}_{1}}$. Поэтому $С_{1}^{4} \leqslant \frac{6}{c}$. Получаем противоречие. Следовательно, верна

Теорема 6. Матрица системы (2.7)–(2.9) в случае представления элементов полей потенциалом двойного слоя имеет ранг не меньше двух, если две точки наблюдения различны.

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

Однако можно привести два примера однозначной разрешимости системы (2.7)–(2.9) при N = 3 и N = 4.

Пример 1: N = 3.

Если число точек наблюдения равно трем, то матрица системы (2.7)–(2.9) для потенциала простого слоя приобретает вид

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

Запишем выражение для детерминанта матрицы (2.21):

Выясним, может ли он быть равен нулю. Имеем

Числа a и b всегда больше нуля при построении S-аппроксимаций. Поэтому последнее равенство может выполняться только в том случае, когда эти числа равны. Но они не могут быть равны, согласно теореме 5. Следовательно, в данном частном случае трехмерной системы уравнений решение определяется однозначно.

Пример 2: N = 4.

1. Рассмотрим систему (2.7)–(2.9) для четырех точек наблюдения при аппроксимации элементов поля с помощью потенциала простого слоя. Матрица системы имеет в этом случае вид

Выберем координаты точек наблюдения так, чтобы элементы матрицы ${{A}_{{44}}}$, обозначенные через ${{d}_{i}},i = 1,...,4,$ были одинаковыми. Почему этого можно добиться?

Рассмотрим поверхность $C = \frac{z}{{{{{({{z}^{2}} + {{r}^{2}})}}^{{3/2}}}}}$ в пространстве координат x, y, z. Выразим ${{r}^{2}}$ через z явно:

Производная функции $f(z)$ положительна при $z < \frac{1}{{{{{(3C{\kern 1pt} ')}}^{{3/4}}}}}$, $C{\kern 1pt} ' = {{C}^{{2/3}}}$.

Функция $f(z)$ является непрерывной функцией своего аргумента и $f(z) \geqslant 0,$ если

Отсюда следует, что при выполнении условий

существуют такие ${{r}_{{1k}}},{{r}_{{2k}}}$, что точки $({{z}_{1}} + {{z}_{k}} + 2H,{{r}_{{1k}}}),$ $({{z}_{2}} + {{z}_{k}} + 2H,{{r}_{{2k}}})$, k = 3, 4, лежат на указанной выше поверхности, и матрица системы имеет вид (2.22). Можно представить выпуклый четырехугольник, в вершинах которого расположены четыре точки наблюдения, z-координаты этих точек различны, но суммы ${{z}_{i}} + {{z}_{k}}$ могут быть одинаковыми при каких-то значениях индексов (но не при всех: четыре точки на плоскости не могут располагаться попарно на одинаковых расстояниях). Вычислим детерминант этой матрицы с помощью формулы Лапласа разложения по двумерным минорам первых двух строк: ${{\varepsilon }_{i}}$ – знак перестановки, соответствующий i-му минору, ${{\bar {M}}_{i}}$ – дополнительный минор.

В итоге получаем

Согласно теореме 5, первое слагаемое всегда положительно, что касается второго слагаемого, то оно тоже положительно в силу соотношения

Таким образом, в данном частном случае четырех точек наблюдения решение системы (2.4) определяется однозначно.

2. Предположим теперь, что четыре точки в трехмерном пространстве расположены таким образом, что матрица системы (2.7)–(2.9) в случае потенциала простого слоя является тёплицевой:

Можно считать, что точки наблюдения находятся в вершинах ромба с углом в 60 градусов. Для определителя (2.21) получаем следующее выражение:

Считаем, что a > b > c (выполнения этого условия можно добиться с помощью указанного выше расположения точек наблюдения). Имеем

Таким образом, мы показали, что в рассмотренном случае матрица системы (2.4) не вырождена, и, следовательно, система имеет единственное решение.

ЗАКЛЮЧЕНИЕ

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

Авторы выражают глубокую благодарность А.С. Леонову за полезные рекомендации и внимание к работе.