On the construction of an optimal network of observation points when solving inverse lin-ear problems of gravimetry and magnetometry

Толық мәтін

О построении оптимальной сети точек наблюдений при решении обратных линейных задач гравиметрии и магнитометрии [1]

ВВЕДЕНИЕ

При интерпретации геофизических данных часто применяются методы, базирующиеся на различных интегральных представлениях магнитного и гравитационного полей [1-6]. Представления полезных сигналов в виде свертки известной функции с подлежащей определению поверхностной, объемной и даже линейной плотностью распределения масс, магнитных диполей и т.п. позволяют редуцировать обратные задачи геофизики по поиску параметров геологической среды к решению систем линейных алгебраических уравнений (СЛАУ) с приближенно заданными матрицей и правыми частями [1]. В наших предыдущих работах [7, 8] были рассмотрены проблемы единственности решения СЛАУ с симметрическими положительно полуопределенными матрицами в локальном и региональном вариантах, а также сформулированы критерии вырожденности СЛАУ, возникающих при использовании метода интегральных уравнений [4, 5].

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

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

$\begin{array}{c}V\left(M\right)={\displaystyle \sum _{l=1}^L{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{{\rho }_{\mathrm{1,}l}({\xi }_1,{\xi }_2)d{\xi }_{1}d{\xi }_2}{\sqrt{{\left(x-{\xi }_1\right)}^2+{\left(y-{\xi }_2\right)}^2+{\left(z-H_l\right)}^2}}+}}}\\ +{\displaystyle \sum _{l=1}^L{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{{\rho }_{\mathrm{2,}l}\left({\xi }_1,{\xi }_2\right)\left(z-H_l\right)d{\xi }_{1}d{\xi }_2}{{\left[\sqrt{{\left(x-{\xi }_1\right)}^2+{\left(y-{\xi }_2\right)}^2+{\left(z-H_l\right)}^2}\right]}^3}}}},\\ M=\left(x,y,z\right),\widehat{x}=\left(x,y\right),\xi =\left({\xi }_1,{\xi }_2,{\xi }_3\right),\widehat{\xi }=\left({\xi }_1,{\xi }_2\right).\end{array}$ (1)

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

$\begin{array}{c}-\frac{\partial V}{\partial z}\left(M\right)={\displaystyle \sum _{l=1}^L{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{{\rho }_{\mathrm{1,}l},\left(\xi \right)\left(z-H_l\right)zd{\xi }_{1}d{\xi }_2}{{\left[\sqrt{{\left(x-{\xi }_1\right)}^2+{\left(y-{\xi }_2\right)}^2+{\left(z-H_l\right)}^2}\right]}^3}+}}}\\ +{\displaystyle \sum _{l=1}^L{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{{\rho }_{\mathrm{2,}l}\left(\xi \right)\left(2{\left(z-H_l\right)}^2-{\left(x-{\xi }_1\right)}^2-{\left(y-{\xi }_2\right)}^2\right)d{\xi }_{1}d{\xi }_2}{{\left[\sqrt{{\left(x-{\xi }_1\right)}^2+{\left(y-{\xi }_2\right)}^2+{\left(z-H_l\right)}^2}\right]}^5}}}},M=\left(x,y,z\right).\end{array}$ (2)

Вторая производная потенциала (например, производные по z и по x) может быть выражена следующим интегралом:

$\begin{array}{c}-\frac{{\partial }^{2}V}{\partial z\partial x}\left(M\right)=3{\displaystyle \sum _{l=1}^L{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{{\rho }_{\mathrm{1,}l},\left(\xi \right)\left(x-{\xi }_1\right)\left(z-H_l\right)d{\xi }_{1}d{\xi }_2}{{\left[\sqrt{{\left(x-{\xi }_1\right)}^2+{\left(y-{\xi }_2\right)}^2+{\left(z-H_l\right)}^2}\right]}^5}+}}}\\ +3{\displaystyle \sum _{l=1}^L{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{{\rho }_{\mathrm{2,}l}\left(\xi \right)\left(x-{\xi }_1\right)\left({\left(x-{\xi }_1\right)}^2+{\left(y-{\xi }_2\right)}^2-4{\left(z-H_l\right)}^2\right)d{\xi }_{1}d{\xi }_2}{{\left[\sqrt{{\left(x-{\xi }_1\right)}^2+{\left(y-{\xi }_2\right)}^2+{\left(z-H_l\right)}^2}\right]}^7}}}},M=\left(x,y,z\right).\end{array}$ (3)

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

$-\frac{\partial V\left(M_i\right)}{\partial z}\equiv f_i={\displaystyle \sum _{l=1}^L{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\left({\rho }_{\mathrm{1,}l}\left(\widehat{\xi }\right)Q_{\mathrm{1,}l}^{(i)}\left(\widehat{\xi }\right)+{\rho }_{\mathrm{2,}l}\left(\widehat{\xi }\right)Q_{\mathrm{2,}l}^{(i)}\left(\widehat{\xi }\right)\right)d\widehat{\xi }}}},i=\mathrm{1,2,...,}N.$ (4)

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

$\Omega \left(\rho \right)={\displaystyle \sum _{l=1}^L{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\left({\rho }_{\mathrm{1,}l}{}^2\left(\widehat{\xi }\right)+{\rho }_{\mathrm{2,}l}{}^2(\left(\widehat{\xi }\right)\right)d\widehat{\xi }}}}=\underset{\rho }{\text{min}},$ (5)

$f_{i,\delta }-{\displaystyle \sum _{l=1}^L{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\left({\rho }_{\mathrm{1,}l}\left(\widehat{\xi }\right)Q_{\mathrm{1,}l}^{(i)}\left(\widehat{\xi }\right)+{\rho }_{\mathrm{2,}l}\left(\widehat{\xi }\right)Q_{\mathrm{2,}l}^{(i)}\left(\widehat{\xi }\right)\right)d\widehat{\xi }}}}=\mathrm{0,}i=\mathrm{1,2,}\dots ,N,$ (6)

находим, что искомые функции должны иметь вид [4, 5, 9]:

$\begin{array}{c}{\rho }_{\mathrm{1,}l}{}^{(a)}\left(\widehat{\xi }\right)={\tilde{\rho }}_{\mathrm{1,}l}\left(\widehat{\xi },\lambda \right),{\rho }_{\mathrm{2,}l}{}^{(a)}\left(\widehat{\xi }\right)={\tilde{\rho }}_{\mathrm{2,}l}\left(\widehat{\xi },\lambda \right),\\ {\tilde{\rho }}_{\mathrm{1,}l}\left(\widehat{\xi },\lambda \right)={\displaystyle \sum _{i=1}^N{\lambda }_i{Q}_{\mathrm{1,}l}^{(i)}\left(\widehat{\xi }\right),{\tilde{\rho }}_{\mathrm{2,}l}\left(\widehat{\xi },\lambda \right)={\displaystyle \sum _{i=1}^N{\lambda }_i{Q}_{\mathrm{2,}l}^{(i)}\left(\widehat{\xi }\right)}},l=\mathrm{1,2,...,}L.\end{array}$ (7)

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

$\mathbf{А}\lambda =f_{\delta },$ (8)

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

$a_{ij}={\displaystyle \sum _{l=1}^L{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\left(Q_{\mathrm{1,}l}^{(i)}\left(\widehat{\xi }\right)Q_{\mathrm{1,}l}^{(j)}\left(\widehat{\xi }\right)+Q_{\mathrm{2,}l}^{(i)}\left(\widehat{\xi }\right)Q_{\mathrm{2,}l}^{(j)}\left(\widehat{\xi }\right)\right)d\widehat{\xi }},}}1\le i\le N,1\le j\le N.$ (9)

Элементы а_ij матрицы А при использовании интегральных представлений (2) и (3) могут быть вычислены явно с помощью интеграла Пуассона. Например, в случае представления вертикальной компоненты гравитационного поля получим выражения:

$\begin{array}{c}{a}_{ij}=2\pi {\displaystyle \sum _{l=1}^L\left\{\frac{z_i+z_j-2H_l}{{\left(\sqrt{{\left(z_i+z_j-2H_l\right)}^2+{\left(x_i-x_j\right)}^2+{\left(y_i-y_j\right)}^2}\right)}^3}\right.}-\\ -\left.\frac{\left(z_i+z_j-2H_l\right)\left(9\left[{\left(x_i-x_j\right)}^2+{\left(y_i-y_j\right)}^2\right]-6{\left(z_i+z_j-2H_l\right)}^2\right)}{{\left(\sqrt{{\left(z_i+z_j-2H_l\right)}^2+{\left(x_i-x_j\right)}^2+{\left(y_i-y_j\right)}^2}\right)}^7}\right\},\\ 1\le i\le N,1\le j\le N.\end{array}$ (10)

По найденным из решения системы (8)–(10) множителям λ_i, i = 1, 2, …, N, можно далее определить высшие производные потенциала, выполнить аналитическое продолжение гравитационного поля и т.д. Для приложений очень важно иметь возможность пересчета значений поля на некоторую регулярную сеть наблюдений. В настоящей работе акцент делается на построении такой сети наблюдений, чтобы матрица системы (8) была невырожденной.

Если какая-либо высшая производная гравитационного потенциала аппроксимируется потенциалом простого слоя, распределенного на одной или нескольких горизонтальных плоскостях, то матрица системы (8)–(10) принимает вид

$\begin{array}{l}A=\left[\begin{array}{ccc}\frac{1}{4{(z_1+H)}^2}& \frac{z_1+z_2+2H}{{\left({\left(z_1+z_2+2H\right)}^2+r_{12}^2\right)}^{3/2}}\mathrm{...}& \frac{z_1+z_N+2H}{{\left({\left(z_1+z_N+2H\right)}^2+r_{1N}^2\right)}^{3/2}}\\ \mathrm{...}& \mathrm{...}& \mathrm{...}\\ \frac{z_1+z_N+2H}{{\left({\left(z_1+z_N+2H\right)}^2+r_{1N}^2\right)}^{3/2}}& \frac{z_2+z_N+2H}{{\left({\left(z_2+z_N+2H\right)}^2+r_{2N}^2\right)}^{3/2}}\mathrm{...}& \frac{1}{4{(z_N+H)}^2}\end{array}\right].\\ \end{array}$ (11)

Матрица системы (11) симметрическая, все ее элементы, как легко видеть, неотрицательные.

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

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

 

Фиг. 1.Проекции областей расположения двух групп точек наблюдения на горизонтальную плоскость.

 

Фиг.2.Схематическое изображение процесса добавления точки в оптимальную сеть.

 

Пример 1. Число точек наблюдений N = 3M. Если число точек наблюдения равно трем, то матрица системы (8)–(10) для потенциала простого слоя приобретает вид

$\begin{array}{c}{A}_{33}=\left[\begin{array}{ccc}a& b& b\\ b& a& b\\ b& b& a\end{array}\right],\\ a=\frac{1}{4{(z_1+H)}^2}=\frac{1}{4{(z_2+H)}^2}=\frac{1}{4{(z_3+H)}^2};\\ b=\frac{z_1+z_2+2H}{{\left({\left(z_1+z_2+2H\right)}^2+r_{12}^2\right)}^{3/2}}=\frac{z_1+z_3+2H}{{\left({\left(z_1+z_3+2H\right)}^2+r_{13}^2\right)}^{3/2}}=\\ =\frac{z_2+z_3+2H}{{\left({\left(z_2+z_3+2H\right)}^2+r_{23}^2\right)}^{3/2}}.\end{array}$ (12)

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

$\frac{z}{{\left(z^2+r^2\right)}^{3/2}}=C=\text{const}>0.$ (13)

 В [8] доказано, что в данном частном случае трехмерной системы уравнений решение определяется однозначно. Если разделить точки наблюдения на группы по три точки в каждой и принять, что z-координаты в каждой такой подсистеме у всех элементов одинаковы, то матрица (12) приобретает вид

$\begin{array}{c}{A}_{3N3N}=\left[\begin{array}{ccc}{A}_{11}& B_{12}& \mathrm{...}{B}_{1N}\\ \mathrm{...}& \mathrm{...}& \mathrm{...}\\ B_{N1}& \mathrm{...}& A_{NN}\end{array}\right],A_{11}=\left[\begin{array}{ccc}a& b& b\\ b& a& b\\ b& b& a\end{array}\right],B_{12}=\left[\begin{array}{ccc}{c}_{11}& c_{12}& c_{13}\\ c_{21}& c_{22}& c_{23}\\ c_{31}& c_{32}& c_{33}\end{array}\right],\\ a=\frac{1}{4{\left(z_1+H\right)}^2}=\frac{1}{4{\left(z_2+H\right)}^2}=\frac{1}{4{\left(z_3+H\right)}^2};\\ b=\frac{z_1+z_2+2H}{{\left({\left(z_1+z_2+2H\right)}^2+r_{12}^2\right)}^{3/2}}=\frac{z_1+z_3+2H}{{\left({\left(z_1+z_3+2H\right)}^2+r_{13}^2\right)}^{3/2}}=C=\\ =\frac{z_2+z_3+2H}{{\left({\left(z_2+z_3+2H\right)}^2+r_{23}^2\right)}^{3/2}};c_{i1}=\frac{z_i+z_4+2H}{{\left({\left(z_i+z_4+2H\right)}^2+r_{i4}^2\right)}^{3/2}},\\ c_{i2}=\frac{z_i+z_5+2H}{{\left({\left(z_i+z_5+2H\right)}^2+r_{i5}^2\right)}^{3/2}},c_{i3}=\frac{z_i+z_6+2H}{{\left({\left(z_i+z_6+2H\right)}^2+r_{i6}^2\right)}^{3/2}},i=\mathrm{1,2,3.}\end{array}$ (14)

В (14) блоки $B_{ij}$ и $A_{ii}$, i, j =1,..,N, имеют вид, аналогичный $A_{11},B_{12}$. Элементы блоков $B_{ij}$ по модулю существенно меньше, чем элементы диагональных блоков, поэтому можно выбрать такое расположение точек наблюдения, что определитель всей матрицы будет отличен от нуля, если отличны от нуля определители диагональных подматриц, поскольку [14] определитель блочно-диагональной матрицы с нулевыми внедиагональными блоками равен произведению определителей диагональных блоков. В нашем случае внедиагональные блоки будут вносить пренебрежимо малый вклад в определитель матрицы всей системы (размерность которой – 3N), а сам по себе определитель является непрерывной функцией элементов матрицы. Можно также вспомнить известное утверждение о блочных матрицах [14]:

$\det \left[\begin{array}{cc}{A}_{}& B_1\\ B_2& C_{}\end{array}\right]=\det A_{}\cdot \det \left[C_{}-B_1{A}_{}^{-1}{B}_2\right].$ (15)

В (15) латинскими буквами обозначены блочные матрицы. Если элементы внедиагональных блоков матрицы (15) много меньше, чем элементы диагональных (по модулю), то определитель второго множителя в (15) мало отличается от определителя блока С (поскольку в произведении трех матриц будут фигурировать элементы, по меньшей мере, второго порядка малости по сравнению с элементами матрицы С). Для более строгого обоснования справедливости нашего утверждения напомним, что определитель суммы двух матриц А и В равен сумме определителей матриц, состоящих из всех возможных комбинаций столбцов матрицы А и матрицы В. Если столбцы матрицы В представляют собой столбцы матрицы А, умноженные на малое число $\epsilon$, то можно написать

$\det \left(A+B\right)=\det A\cdot {\displaystyle \sum _{k=0}^N{C}_N^k}{\epsilon }^k=\det A\cdot {\left(1+\epsilon \right)}^N,C_N^k=\frac{N!}{k!\left(N-k\right)!}.$

Таким образом, определитель матрицы А+В будет не слишком сильно отличаться от определителя А, если $\epsilon$ достаточно мало.

Пример 2. Число точек наблюдений N = 4M. Рассмотрим сначала, как и в работе [8], систему (8)–(10) для четырех точек наблюдения при аппроксимации элементов поля с помощью потенциала простого слоя. Матрица системы имеет в этом случае вид

$\begin{array}{c}{A}_{44}=\left[\begin{array}{cccc}a& c& d_1& d_2\\ c& b& d_3& d_4\\ d_1& d_3& a_1& c_1\\ d_2& d_4& c_1& b_1\end{array}\right],\\ a=\frac{1}{4{(z_1+H)}^2};c=\frac{z_1+z_2+2H}{{\left({\left(z_1+z_2+2H\right)}^2+r_{12}^2\right)}^{3/2}};\\ b=\frac{1}{4{(z_2+H)}^2};a_1=\frac{1}{4{(z_3+H)}^2};\\ b_1=\frac{1}{4{(z_4+H)}^2};c_1=\frac{z_3+z_4+2H}{{\left({\left(z_3+z_4+2H\right)}^2+r_{34}^2\right)}^{3/2}};\\ d_1=\frac{z_1+z_3+2H}{{\left({\left(z_1+z_3+2H\right)}^2+r_{13}^2\right)}^{3/2}};d_2=\frac{z_1+z_4+2H}{{\left({\left(z_1+z_4+2H\right)}^2+r_{14}^2\right)}^{3/2}};\\ d_3=\frac{z_2+z_3+2H}{{\left({\left(z_2+z_3+2H\right)}^2+r_{23}^2\right)}^{3/2}};d_4=\frac{z_2+z_4+2H}{{\left({\left(z_4+z_4+2H\right)}^2+r_{24}^2\right)}^{3/2}}.\end{array}$ (16)

Выберем координаты точек наблюдения так, чтобы элементы матрицы $A_{44}$, обозначенные через $d_i,i=\mathrm{1,...,4,}$ были одинаковыми. Ранее мы показали, как этого можно добиться, принимая во внимание свойства функции, связывающей z-координаты точек наблюдений с расстоянием между проекциями этих точек на горизонтальную плоскость (см. [8]). Решение системы (8)$–$(10) в случае 4 точек наблюдения, расположенных на одной горизонтальной плоскости в вершинах ромба, определяется однозначно. Так же, как и в примере 1, перейдем к блочно-диагональной матрице для системы 4M точек наблюдения. Покажем, как можно осуществить подобный переход путем последовательного увеличения числа пунктов измерений.

Предположим, что нам даны 8 точек наблюдений, и матрица системы имеет вид

$\begin{array}{c}{A}_{88}=\left[\begin{array}{cc}{A}_{44}& B_{44}\\ B_{44}& C_{44}\end{array}\right],A_{44}=\left[\begin{array}{cccc}a& c& d_1& d_2\\ c& b& d_3& d_4\\ d_1& d_3& a_1& c_1\\ d_2& d_4& c_1& b_1\end{array}\right],\\ a=\frac{1}{4{\left(z_1+H\right)}^2};c=\frac{z_1+z_2+2H}{{\left({\left(z_1+z_2+2H\right)}^2+r_{12}^2\right)}^{3/2}};\\ b=\frac{1}{4{\left(z_2+H\right)}^2};a_1=\frac{1}{4{\left(z_3+H\right)}^2};\\ b_1=\frac{1}{4{\left(z_4+H\right)}^2};c_1=\frac{z_3+z_4+2H}{{\left({\left(z_3+z_4+2H\right)}^2+r_{34}^2\right)}^{3/2}};\\ d_1=\frac{z_1+z_3+2H}{{\left({\left(z_1+z_3+2H\right)}^2+r_{13}^2\right)}^{3/2}};d_2=\frac{z_1+z_4+2H}{{\left({\left(z_1+z_4+2H\right)}^2+r_{14}^2\right)}^{3/2}};\\ d_3=\frac{z_2+z_3+2H}{{\left({\left(z_2+z_3+2H\right)}^2+r_{23}^2\right)}^{3/2}};d_4=\frac{z_2+z_4+2H}{{\left({\left(z_4+z_4+2H\right)}^2+r_{24}^2\right)}^{3/2}}.\end{array}$ (17)

Матрицы $B_{44}$ и $C_{44}$ – такие 4х4 – матрицы, что определитель всей матрицы $A_{88}$ отличен от нуля. Этого можно добиться следующим образом.

Сначала выбираем такие четыре точки, чтобы определитель четырехмерной матрицы $A_{44}$ не был равен нулю, а затем учитываем справедливость выражения [14]:

$\det \left[\begin{array}{cc}{A}_{44}& B_{44}\\ B_{44}& C_{44}\end{array}\right]=\det A_{44}\cdot \det \left[C_{44}-B_{44}{A}_{44}^{-1}{B}_{44}\right].$

Координаты четырех точек «второй очереди» можно задать так, чтобы расстояния между ними и точками «первой очереди» были существенно больше, чем расстояния между точками, принадлежащими каждой группе в отдельности. Тогда определитель второго множителя будет также отличен от нуля. Далее действуем так же, как и в примере 1.

Таким образом, нами доказана

Теорема 1. В случае N = 3M и N = 4M (N – число точек наблюдений, M – натуральное число) существуют такие конфигурации сети наблюдений, что матрица системы (8)–(10) является невырожденной.

Пример 3. Случай произвольного N. Допустим теперь, что точки наблюдения расположены таким образом, что выполняется достаточное условие невырожденности матрицы [14]:

$\left|d_{ii}\right|>{\displaystyle \sum _{j\ne i}^N\left|d_{ij}\right|},i=\mathrm{1,2,...}N\Rightarrow \det D\ne 0.$ (18)

Слева в (18) стоят модули диагональных элементов некоторой матрицы, а справа $–$ сумма модулей элементов каждой строки в отдельности.

Рассмотрим диагональные и внедиагональные элементы матрицы (8):

$a=\frac{1}{4{(z_i+H)}^2};c=\frac{z_i+z_j+2H}{{\left({\left(z_i+z_j+2H\right)}^2+r_{ij}^2\right)}^{3/2}},i,j=\mathrm{1,...,}N.$

Если вертикальные координаты всех точек одинаковы, то достаточное условие невырожденности матрицы (18) будет выполнено в том случае, когда

$r_{ij}^2\ge \left({\left(N-1\right)}^{2/3}-1\right)\cdot {\left(z+H\right)}^2,z_i\equiv z;i,j=\mathrm{1,...,}N.$ (19)

Итак, верна

Теорема 2. При выполнении условия (19) решение системы (8)–(10) определяется однозначно.

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

Пример 4. Поверхность уровня гармонической функции. Рассмотрим систему (8)$–$(10) для произвольного числа точек наблюдений М и предположим, что точки наблюдения находятся на поверхности уровня гармонической функции, равной в каждой точке вертикальной компоненте гравитационного потенциала:

$С=\frac{z_i+z_j+2H}{{\left({\left(z_i+z_j+2H\right)}^2+r_{ij}^2\right)}^{3/2}},i,j=\mathrm{1,...,}N.$ (20)

Исследуем вопрос однозначной разрешимости системы (8)$–$(10) в этом случае.

Если z-координаты точек наблюдения известны, то из (20) следует, что расстояния между проекциями точек на плоскости можно найти по следующей формуле:

$\begin{array}{c}{r}_{ij}^2={\left(\frac{z_i+z_j+2H}{Ñ}\right)}^{2/3}-{\left(z_i+z_j+2H\right)}^2,\\ i,j=\mathrm{1,...,}N.\end{array}$ (21)

В (21) предполагается, что поле создается только одним простым слоем, расположенным на глубине Н.

Пусть матрица СЛАУ имеет следующий вид:

$A_{NN}=\left[\begin{array}{l}\begin{array}{ccc}{a}_{11}& 1& 1\\ 1& a_{22}& 1\\ 1& 1& \mathrm{...}\end{array}\begin{array}{ccc}\mathrm{...}& 1& 1\\ 1& \mathrm{...}& 1\\ 1& 1& 1\end{array}\\ \begin{array}{ccc}\mathrm{...}& \mathrm{...}& \mathrm{...}\\ 1& 1& \mathrm{...}\\ 1& 1& \mathrm{...}\end{array}\begin{array}{ccc}\mathrm{...}& \mathrm{...}& 1\\ \mathrm{...}& a_{N-1N-1}& 1\\ \mathrm{...}& \mathrm{...}& a_{NN}\end{array}\end{array}\right].$ (22)

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

Теорема 3. Определитель матрицы равен

$\det A_{NN}=\left(a_{11}+{\displaystyle \sum _{j=2}^N\frac{a_{11}-1}{a_{jj}-1}}\right)\times \left(a_{22}-1\right)\cdot \left(a_{33}-1\right)\cdot \left(a_{NN}-1\right).$ (23)

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

$\det A_{NN}=\left|\begin{array}{l}\begin{array}{ccc}{a}_{11}-{\displaystyle \sum _{j=2}^N\frac{1-a_{11}}{a_{jj}-1}}& 1-a_{11}& \mathrm{...}1-a_{11}\\ 0& a_{22}-1& \mathrm{...}0\\ 0& 0& a_{N-1N-1}0\end{array}\\ 00\mathrm{...}{a}_{NN}\end{array}\right|.$

Перемножение диагональных элементов приводит нас к утверждению теоремы.

Доказанная нами теорема 3 может быть полезна при построении сети точек наблюдений, которая гарантировала бы невырожденность соответствующей СЛАУ при интерпретации гравитационных и магнитных данных. А именно, следует смоделировать такую сеть, чтобы точки наблюдения находились на поверхности уровня некоторой гармонической функции, значениями которой в этих точках являются элементы матрицы СЛАУ. Однако не все так просто, как кажется на первый взгляд. Проблема построения сети заключается в том, что элементы матрицы $–$ это не совсем абстрактные вещественные числа. Элементы СЛАУ содержат в себе информацию о расстояниях между точками сети наблюдений как в трехмерном пространстве, так и на плоскости (тогда рассматриваются проекции соответствующих точек). Поэтому «в лоб» применять формулу (21) для нахождения расстояний между проекциями точек на плоскости, если известны z-координаты этих точек, нельзя: необходимо сформулировать некоторую вариационную задачу, решение которой позволило бы достичь минимального отклонения координат реальных точек от расчетных. Кроме того, диагональные элементы матрицы должны быть все либо больше внедиагональных, либо меньше. Если сумма в первом сомножителе (23) знакопеременная, то определитель может при некоторых значениях диагональных элементов обращаться в нуль. Для того чтобы избежать такой ситуации, следует строить оптимальную сеть наблюдений таким образом, чтобы определитель в (23) был гарантированно отличен от нуля.

Пример 5. Случай кососимметрических матриц. В тех случаях, когда необходимо построить математические модели высших производных гравитационного или магнитного потенциала (см. формулу (3)), матрица системы линейных алгебраических уравнений оказывается кососимметрической. На диагонали стоят нули, а внедиагональные элементы «при отражении от главной диагонали» меняют знак. Если элементы матрицы (8) продифференцировать по $x_i$, например, то мы как раз и получим описанную выше систему линейных алгебраических уравнений с кососимметрической матрицей.

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

Докажем следующую теорему.

Теорема 4. Определитель матрицы

$A_{NN}=\left[\begin{array}{l}\begin{array}{ccc}0& -1& -1\\ 1& 0& -1\\ 1& 1& \mathrm{...}\end{array}\begin{array}{ccc}\mathrm{...}& -1& -1\\ -1& \mathrm{...}& -1\\ -1& -1& -1\end{array}\\ \begin{array}{ccc}\mathrm{...}& \mathrm{...}& \mathrm{...}\\ 1& 1& \mathrm{...}\\ 1& 1& \mathrm{...}\end{array}\begin{array}{ccc}\mathrm{...}& \mathrm{...}& -1\\ \mathrm{...}& 0& -1\\ \mathrm{...}& \mathrm{...}& 0\end{array}\end{array}\right]$ (24)

равен 1 в случае четного N и нулю в противном случае.

Доказательство. Вычтем первый столбец из последующих, а затем $–$ последнюю строку из всех остальных. Получим матрицу

$A_{NN}=\left[\begin{array}{l}\begin{array}{ccc}0& -2& -2\\ 0& -1& -2\\ 0& 0& \mathrm{...}\end{array}\begin{array}{ccc}\mathrm{...}& -2& -1\\ -2& \mathrm{...}& -1\\ -2& -2& -1\end{array}\\ \begin{array}{ccc}\mathrm{...}& \mathrm{...}& \mathrm{...}\\ 0& 0& \mathrm{...}\\ 1& 0& \mathrm{...}\end{array}\begin{array}{ccc}\mathrm{...}& \mathrm{...}& \mathrm{...}-1\\ \mathrm{...}& -1& -1\\ \mathrm{...}& \mathrm{...}& -1\end{array}\end{array}\right].$ (25)

Разложим определитель матрицы (25) по первому столбцу.

Алгебраическое дополнение последнего элемента этого столбца равно:

${\Delta }_{N1}={\left(-1\right)}^{N-1}.$

Поэтому определитель всей матрицы равен $\det A_{NN}={\left(-1\right)}^{N-1+N+1}=1.$

Если число точек наблюдений четное, то определитель равен 1, в противном случае он равен нулю.

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

$C=\frac{(x_i-x_j)\left(z_i+z_j+2H\right)}{{\left[\sqrt{{\left(x_i-x_j\right)}^2+{\left(y_i-y_j\right)}^2+{\left(z_i+z_j+2H\right)}^2}\right]}^5},i,j=\mathrm{1,...}N.$

Как и ранее, Н – это глубина залегания эквивалентных источников, распределенных на плоскости (простого слоя, см. (2) и (3)).

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

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

$\begin{array}{l}A=\left[\begin{array}{ccc}\frac{6}{16{(z_1+H)}^4}& \mathrm{...}& \frac{\left(z_1+z_N+2H\right)\cdot \left(6{\left(z_1+z_N+2H\right)}^2-9r_{1N}^2\right)}{{\left({\left(z_1+z_N+2H\right)}^2+r_{1N}^2\right)}^{7/2}}\\ \mathrm{...}& \mathrm{...}& \mathrm{...}\\ \frac{\left(z_1+z_N+2H\right)\cdot \left(6{\left(z_1+z_N+2H\right)}^2-9r_{1N}^2\right)}{{\left({\left(z_1+z_N+2H\right)}^2+r_{1N}^2\right)}^{7/2}}& \mathrm{...}& \frac{6}{16{(z_N+H)}^4}\end{array}\right].\\ \end{array}$ (26)

Если применить теоремы 2 и 3, то мы получим невырожденную матрицу СЛАУ. Теорема 2 утверждает, что матрица (26) будет иметь отличный от нуля определитель, если выполняется следующее условие:

$r_{ij}^2\ge \left({\left(N-1\right)}^{2/7}-1\right)\cdot {\left(z+H\right)}^2,z_i\equiv z;i,j=\mathrm{1,...,}N.$ (27)

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

$\frac{\left(z_i+z_j+2H\right)\cdot \left(6{\left(z_i+z_j+2H\right)}^2-9r_{ij}^2\right)}{{\left({\left(z_i+z_j+2H\right)}^2+r_{ij}^2\right)}^{7/2}}=С=\text{const}.$ (28)

В (28) Н $–$ глубина залегания двойного слоя.

2. АЛГОРИТМ ПОСТРОЕНИЯ ОПТИМАЛЬНОЙ СЕТИ ТОЧЕК

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

При «конструировании» сети наблюдений можно пойти двумя путями.

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

$\begin{array}{c}\Phi \left[Z\right]\equiv {{\displaystyle \sum _{i,j}^N\left(r_{ij}^2-f_{ij}(Z)\right)}}^2,\\ f_{ij}={\left(\frac{(z_i+z_j+2H)}{C}\right)}^{2/3}-{\left(z_i+z_j+2H\right)}^2,i,j=\mathrm{1...}N.\end{array}$ (29)

Поставим вариационную задачу: найти $Z^{*}$

$\begin{array}{l}\underset{\text{}Z}{Z^{*}:\min }\Phi \left[Z\right]+\alpha \Omega \left[Z\right],Z={\left(z_1,z_2\mathrm{,...,}{z}_N\right)}^T,\Phi \left[Z\right]={{\displaystyle \sum _{i,j}^N\left(r_{ij}^2-f_{ij}(Z)\right)}}^2,\\ f_{ij}={\left(\frac{(z_i+z_j+2H)}{C}\right)}^{2/3}-{\left(z_i+z_j+2H\right)}^2,i,j=\mathrm{1...}N,\Omega \left[Z\right]={\displaystyle \sum _{i=1}^N{z}_i^2}.\end{array}$ (30)

Минимум в задаче (30) можно искать одним из градиентных методов [2, 3].

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

$z_1=z_2=z_3,C=\frac{z_1+z_2+2H}{{\left(\sqrt{{\left(z_1+z_2+2H\right)}^2+r_{12}^2}\right)}^3}.$ (31)

Далее будем добавлять к выбранным нами трем точкам последовательно по одному пункту наблюдений, учитывая, что каждый новый пункт должен находиться либо на поверхности уровня (31), либо вблизи нее. Четвертая точка не может лежать в той же горизонтальной плоскости, как это уже отмечалось ранее. Поэтому выберем $z_4$ так, чтобы выполнялись соотношения (см. теорему 3) и рассмотрим треугольник $A_2{A}_3{A}_4.$:

$z_4<z_1\Rightarrow \frac{1}{z_4^2}>\frac{1}{z_1^2}>C.$ 

Из рисунка легко установить, что выполняются следующие соотношения для расстояний между точками треугольника $A_2{A}_3{A}_4$:

$\begin{array}{c}{y}^2=r_{34}^2-x^2=r_{24}^2-{\left(r_{23}-x\right)}^2\Rightarrow x=\frac{r_{24}^2-r_{34}^2-r_{23}^2}{2r_{23}},\\ y^2=r_{34}^2-{\left(\frac{r_{24}^2-r_{34}^2-r_{23}^2}{2r_{23}}\right)}^2,\overrightarrow{A_{3}O}=\frac{\left(x_2-x_3,y_2-y_3\right)}{r_{23}}\cdot x,\\ \overrightarrow{r_4}=\overrightarrow{r_3}+\overrightarrow{A_{3}O}+y\cdot \frac{\left(y_3-y_2,x_3-x_2\right)}{r_{23}}.\end{array}$ (32)

В (32) через $r_i,i=\mathrm{2,3,4,}$ обозначены радиус-векторы соответствующих точек в некоторой системе координат. Если мы задаем расстояния от четвертой точки до второй и третьей, то из (32) становится ясно, что четвертая точка может находиться только в двух положениях на плоскости, определяемых соотношениями (32) (по обе стороны от отрезка $A_2{A}_3$ ). Таким образом, расположение точки $A_4$ относительно $A_1$ не может быть произвольным и не обязательно должно быть справедливо равенство

$r_{14}^2={\left(\frac{z_1+z_4+2H}{Ñ}\right)}^{2/3}-{\left(z_1+z_4+2H\right)}^2.$ (33)

Мы приходим к необходимости такого подбора параметров С и $z_4$, чтобы расстояние между первой и четвертой точками, определяемое по (33), меньше всего отличалось от расстояния между этими точками, которое вычисляется с учетом фактического расположения четвертой точки из соотношений (32). Все остальные точки сети можно добавлять аналогичным образом, задавая их расстояния на плоскости (фактически, расстояния между их проекциями на исходную плоскость, которой принадлежат первые три точки) до двух предыдущих с номерами $i-\mathrm{1,}i$: $r_{i-\mathrm{1,}i+1},r_{i,i+1}$. У нас на каждом таком шаге будет получаться треугольник $A_{i-1}{A}_i{A}_{i+1}$ и нам нужно будет выбирать $C,z_4\mathrm{,...,}{z}_{i+1}$ так, чтобы соотношения (32) и (33) для всех пар точек не противоречили друг другу.

ЗАКЛЮЧЕНИЕ

Метод S-аппроксимаций позволяет построить за сравнительно короткое время аналитические аппроксимации элементов аномальных потенциальных полей в «первом приближении» по точности. Такая аппроксимация позволяет составить общее представление о природе источников, характере аномалии и уточнить геологическое строение Земли и планет Солнечной системы. Приведенные в настоящей статье теоремы позволяют проектировать оптимальные сети наблюдений, гарантирующие однозначную разрешимость СЛАУ, к которым редуцируются линейные обратные задачи геофизики, а также разрабатывать эффективные алгоритмы отбора спутниковых данных, которые с большой степенью надежности могут применяться при интерпретации разноточной и разнородной информации об источниках физических полей. Если выборки спутниковых данных не удовлетворяют критериям оптимальности сети наблюдений, то мы получаем плохо обусловленные СЛАУ большой и сверхбольшой размерности, и приближенные решения таких систем не могут использоваться для локализации источников аномалий и нахождения пространственного распределения различных элементов полей.

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

 

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

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

I. Stepanova

Moscow State University

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

D. Lukyanenko

Schmidt Institute of Physics of the Earth of the Russian Academy of Sciences

Email: tet@ifz.ru
Ресей, ul. Bolshaya Gruzinskaya, 10, build. 1, Moscow, 123995

I. Kolotov

Schmidt Institute of Physics of the Earth of the Russian Academy of Sciences

Email: tet@ifz.ru
Ресей, ul. Bolshaya Gruzinskaya, 10, build. 1, Moscow, 123995

A. Shchepetilov

Schmidt Institute of Physics of the Earth of the Russian Academy of Sciences

Email: tet@ifz.ru
Ресей, ul. Bolshaya Gruzinskaya, 10, build. 1, Moscow, 123995

A. Yagola

Schmidt Institute of Physics of the Earth of the Russian Academy of Sciences

Email: tet@ifz.ru
Ресей, ul. Bolshaya Gruzinskaya, 10, build. 1, Moscow, 123995

A. Levashov

Schmidt Institute of Physics of the Earth of the Russian Academy of Sciences

Email: tet@ifz.ru
Ресей, ul. Bolshaya Gruzinskaya, 10, build. 1, Moscow, 123995

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

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