Approximation approach in various modifications of the method of linear integral representations

Full Text

1. ВВЕДЕНИЕ

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

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

1. Возрастает объем данных, подлежащих интерпретации.

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

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

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

5. Появляются новые инженерно-эксплуа­тационные, экономические и экологические проблемы.

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

Владимир Николаевич Страхов неоднократно подчеркивал важность применения аппроксимационного подхода к решению любых задач геофизики. Для лучшего понимания сущности метода интегральных представлений, перспектив его развития и совершенствования целесообразно вспомнить основные идеи аппроксимационного подхода, изложенные В.Н. Страховым в его работе «Геофизика и математика» [Страхов, 1999].

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

Fk=Fk(u), k=1, … K,

где u=u(x) – некоторый элемент потенциального поля (потенциал либо какая-то производная от потенциала по координатам), т. е. функция, определенная на поверхности и во внешности Земли и являющаяся функцией из некоторого линейного нормированного пространства Х (обычно Х = Н – гильбертову пространству), а Fk (·) – функционалы над пространством Х.

В случае задач класса I сначала по заданным из наблюдений величинам:

$f_{i,\delta }=f_i+\delta f_i\text{, }i=1,\text{ 2, }\dots \text{, }J$

где 

$f_i={\Theta }_i\left(u\right)$, 

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

$u_{appr}\left(x\right)={\tilde{u}}_a\left(x\right)$, 

т.  е.  аппроксимации  функций  u(x),  желательно  глобальные,  а  затем  используется  фундаментальное  соотношение:

$F_k=F_k\left(u\right)\approx F_k^{\left(a\right)}\text{, }k=1,\text{ 2, }\dots \text{, }K$, 

где $F_k^{\left(a\right)}(•)$ суть  аппроксимации  функционалов $F_k^{}(•)$.  Во  многих  задачах  аппроксимация  функционалов $F_k^{}(•)$ не  требуется,  т.  е.  используются  соотношения:

$\left.F_k\right|\approx F_k\left({\tilde{u}}_a\right).$

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

${\Phi }_m=F_m\left(\rho \right)$

где: $\rho =\rho \left(\xi \right)$ есть  либо  скалярная  функция  (случай  гравитационного  поля),  либо  вектор-функция  (случай  магнитного  поля),  описывающая  распределение  источников  аномального  поля, $\xi$–  координаты  точек  из  области  V  ниже  поверхности  Земли.  Функция  (вектор-функция) $\rho \left(\xi \right)$ считается  принадлежащей  к  функциональному  пространству  Y(V)  (обычно  –  гильбертову),  а ${\Phi }_m(•)$  –  функционалы  над  пространством  Y(V).

В  задачах  класса  II  сначала  по  заданным  величинам $f_{i,\delta }=f_i+\delta f_i,f_i={\Psi }_i\left(\rho \right)$,  строятся  функции:

${\rho }_{appr}\left(\xi \right)={\tilde{\rho }}_a\left(\xi \right)$,

аппроксимирующие  функции $\rho \left(\xi \right)$ в  априорно  принятом  классе  функций  из  Y(V),  а  затем  используется  фундаментальное  соотношение

${\Phi }_m={\Phi }_m\left(\rho \right)\approx {\Phi }_m^{\left(a\right)}\left({\tilde{\rho }}_a\right)$, 

где ${\Phi }_m^{\left(a\right)}(•)$ суть  аппроксимации  функционалов  ${\Phi }_m^{}(•)$.  Во  многих  задачах  аппроксимация  функ­ционалов ${\Phi }_m^{}(•)$ не  требуется,  т.  е.  используется  соотношение:

${\Phi }_m\approx {\Phi }_m\left({\tilde{\rho }}_a\right).$

Ясно,  далее,  что  многие  задачи  гравиметрии  и  магнитометрии  могут  трактоваться  (и,  соответственно,  решаться)  и  как  задачи  класса  I,  и  как  задачи  класса  II. Например,  таковыми  являются  задачи  нахождения  пространственного  распределения  полей.  Важно  подчеркнуть,  что  аппроксимационный  подход,  по  существу,  является  обобщением  и  развитием  классического  метода  подбора  при  нахождении  аппроксимаций  распределений  источников  аномальных  полей  (или,  как  принято  говорить  в  настоящее  время,  –  при  решении  «обратных  задач»).

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

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

Две  версии  метода  S-аппроксимаций  (локальная  и  региональная),  F-  и  R-аппроксимации  являются  вариантами  метода  линейных  интегральных  представлений;  основные  характеристики  этого  метода  изложены  в  более  ранних  работах  авторов  (см.,  напр.  [Страхов,  2000;  2002а;  2002б;  Степанова,  2009а;  2009б;  2008;  2007;  Раевский,  2015а;,  2015б]).  Разработка  новых  способов  решения  обратных  задач  геофизики  является  актуальной  задачей  ввиду  сложности  и  некорректности  их  постановок  (см.,  напр.,  работы  [Глазнев  и  др.,  1989;  Долгаль,  1994;  Кризский,  2011;  Мартышко,  1997]).

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

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

В  методе  F-аппроксимаций  элементы  аномальных  потенциальных  полей  представляются  интегралом  Фурье,  а  R-аппроксимации  получаются  при  так  называемом  лучевом  преобразо­вании.

Методы  F-,  R-  и  S-аппроксимаций  позволяют  получить  решение,  с  помощью  которого  можно  эффективно  строить  линейные  трансформанты  поля,  а  также  использовать  его  в  качестве  нулевого  приближения  для  решения  нелинейной  обратной  задачи  по  локализации  источников.

В  данной  работе  представлены  результаты  построения  аналитических  аппроксимаций  аномального  гравитационного  поля  и  рельефа  в  локальном  варианте  с  использованием  системы  прямоугольных  декартовых  координат  в  трехмерном  и  четырехмерном  пространствах.  При  этом  существенным  является  тот  факт,  что  S-,  F-,  R-аппроксимации  тесно  связаны  друг  с  другом.  Возникает  идея  применить  все  три  описанных  ранее  в  работах  авторов  варианта  метода  линейных  интегральных  представлений  совместно  с  целью  построения  более  качественных  моделей  геоида  и  аномальных  геопотенциальных  полей.

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

Для  того  чтобы  понять,  почему  существует  тесная  связь  между  различными  вариантами  метода  линейных  интегральных  представлений,  необходимо  напомнить,  как  именно  строятся  R-  и  F-аппроксимации.

Как  было  подробно  описано  в  работе  [Степанова,  2009а],  для  функции $f\left(x\right)\in S\left(R^n\right)$,  где $S\left(R^n\right)$ –  пространство  быстро  убывающих  на  бесконечности  непрерывно  дифференцируемых  функций  (точнее  говоря,  для  непрерывно  дифференцируемых  функций,  имеющих  порядок  убывания $O{(1+\sum _{i=1}^n{x}_i^2)}^{-1}$)  –  пространство  Шварца  –  существует  преобразование  Радона:

 $\hat{f}(\omega ,p)=\underset{(\omega ,x)=p}{\overset{}{\int }}f(x\left)dm\right(x),$ (1) 

где: $\omega$ –  единичный  вектор,  dm(x)  –  мера  на  прямой  $(\omega ,x)=p$.

В  двумерном  случае  формула  (1)  принимает  вид:

$\begin{array}{c}\hat{f}(\omega ,p)=\underset{-\infty }{\overset{\infty }{\int }}f(-t\sin s+x_1\cos \phi ,t\cos s+x_2\sin \phi )ds,\\ \omega =(\cos \phi ,\sin \phi ),x=(x_1,x_2).\end{array}$ (2) 

Запишем  основную  формулу  теории  гармонических  функций  для  полупространства,  ограниченного  плоскостью $x_3=0$ (далее  упоминаемой  как  плоскость  «П»)  [Кошляков  и  др.,  1962]:

$\begin{array}{c}V\left(M\right)=\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}\frac{{\rho }_1({\xi }_1,{\xi }_2)d{\xi }_{1}d{\xi }_2}{\sqrt{{(x_1-{\xi }_1)}^2+{(x_2-{\xi }_2)}^2+{x_3}^2}}+\\ +\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}\frac{{\rho }_2({\xi }_1,{\xi }_2)x_{3}d{\xi }_{1}d{\xi }_2}{{\left[\sqrt{{(x_1-{\xi }_1)}^2+{(x_2-{\xi }_2)}^2+{x_3}^2}\right]}^3},\\ M=(x_1,x_2,x_3),\xi =({\xi }_1,{\xi }_2,{\xi }_3).\end{array}$ (3) 

Мы  выбрали  систему  координат  так,  чтобы  плоскость  простого  и  двойного  слоев  задавалась  уравнением ${х}_3=0$.  Тогда  производная  по  х₃  потенциала  V,  взятая  с  обратным  знаком,  будет  иметь  вид:

$\begin{array}{c}-\frac{\partial V}{\partial x_3}\left(M\right)=\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}\frac{{\rho }_1\left(\hat{\xi }\right)x_{3}d\hat{\xi }}{{\left[{(x_1-{\xi }_1)}^2+{(x_2-{\xi }_2)}^2+x_3^2\right]}^{3/2}}+\\ +\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}\frac{{\rho }_2\left(\hat{\xi }\right){(2{x_3}^2-{(x_1-{\xi }_1)}^2-{(x_2-{\xi }_2)}^2)}^{2}d\hat{\xi }}{{\left[{(x_1-{\xi }_1)}^2+{(x_2-{\xi }_2)}^2+x_3^2\right]}^{5/2}},\\ M=(x_1,x_2,x_3),\text{\hspace{0.33em}}\hat{\xi }=({\xi }_1,{\xi }_2).\end{array}$ (4) 

Функции  ${\rho }_1,{\rho }_2$  неизвестны.  Пусть  компоненты  поля  заданы  в  конечном  множестве  точек $M_i=(x_1^{\left(i\right)},x_2^{\left(i\right)},x_3^{\left(i\right)}),i=1,2,...,N$.  Обозначим  подынтегральную  функцию  в  первом  слагаемом  в  (4)  в  точке  $M_i$ через $Q_1^{\left(i\right)},$ а  во  втором  слагаемом  –  через $Q_2^{\left(i\right)}$.  Тогда  получим:

$\begin{array}{c}-\frac{\partial V\left(M_i\right)}{\partial x_3}\equiv f_i=\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}\left({\rho }_1\right(\hat{\xi }\left)Q_1^{\left(i\right)}\right(\hat{\xi })+\\ +{\rho }_2\left(\hat{\xi }\right)Q_2^{\left(i\right)}\left(\hat{\xi }\right))d\hat{\xi },i=1,2,...,N.\end{array}$ (5)

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

Применив  к  обеим  частям  равенства  (5)  преобразование  Радона,  получим:

$\begin{array}{c}{\hat{V}}_{x_3}(\omega ,p)=\underset{-\infty }{\overset{+\infty }{\int }}\left[{\hat{\rho }}_1\right(\omega ,q\left){\hat{Q}}_1^{\left(i\right)}\right(\omega ,p-q)+\\ +{\hat{\rho }}_2(\omega ,q){\hat{Q}}_2^{\left(i\right)}(\omega ,p-q)]dq,\end{array}$ (6)

т.  к.  преобразование  Радона  свертки  функций  равно  свертке  от  произведения  преобразований  Радона  соответствующих  функций.  Величины  (точнее  говоря,  функции) ${\hat{Q}}_1^{\left(i\right)},{\hat{Q}}_2^{\left(i\right)}$ мы  можем  найти:

$\begin{array}{c}\underset{-\infty }{\overset{\infty }{\int }}\frac{x_3}{{\left[\sqrt{{(x_1+t\sin \phi -p\cos \phi )}^2+{(x_2-t\cos \phi -p\sin \phi )}^2+x_3^2}\right]}^3}dt=\\ =\frac{2x_3^{}}{x_3^2+p^2-2p{x}_1\cos \phi -2p{x}_2\sin \phi +{{(x_1\cos \phi +x_2\sin \phi )}^2}_{{}^{}}},\omega =(\cos \phi ,\sin \phi ),\\ \underset{-\infty }{\overset{\infty }{\int }}\frac{(2{x_3}^2-({(x_1+t\sin \phi -p\cos \phi )}^2+{(x_2-t\cos \phi -p\sin \phi )}^2)}{{\left[\sqrt{{(x_1+t\sin \phi -p\cos \phi )}^2+{(x_2-t\cos \phi -p\sin \phi )}^2+x_3^2}\right]}^5}\text{\hspace{0.17em}}dt=\\ =\frac{\partial }{\partial x_3}\left(\frac{2x_3^{}}{x_3^2+p^2-2p{x}_1\cos \phi -2p{x}_2\sin \phi +{{(x_1\cos \phi +x_2\sin \phi )}^2}_{{}^{}}}\right),\omega =(\cos \phi ,\sin \phi ),\end{array}$ (7) 

где ${\omega }_1\xi +{\omega }_2\eta =p$ – прямая, по которой производится интегрирование.

На  практике  компоненты  поля  бывают  заданы  с  некоторой  погрешностью,  поэтому  входной  информацией  являются  значения $f_i,\text{\hspace{0.33em}}\delta$.  С  помощью  решения  вариационной  задачи  (в  общем  виде  вариационная  постановка  описывается  в  статье  [Страхов  и  др.,  1997]):

$\begin{array}{c}\Omega \left(\rho \right)=\underset{0}{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}dq\underset{0}{\overset{2\pi }{\int }}\left({{\hat{\rho }}_1}^2\text{}\right(\omega ,q)+{{\hat{\rho }}_2}^2(\omega ,q\left)\right)dpd\phi =\underset{\rho }{\text{min}},\\ f_{i,\delta }=-\frac{1}{2\pi }\underset{0}{\overset{2\pi }{\int }}d\phi \underset{-\infty }{\overset{\infty }{\int }}\left\{\left[\underset{-\infty }{\overset{\infty }{\int }},\frac{{\hat{Q}}_{1p}^{\left(i\right)}(\omega ,p-q)dp}{p}\right]{\rho }_1(\omega ,q)+\right.\\ +\left.\left[\underset{-\infty }{\overset{\infty }{\int }},\frac{{\hat{Q}}_{2p}^{\left(i\right)}(\omega ,p-q)dp}{p}\right]{\rho }_2(\omega ,q)\right\}dq\end{array}$ (8) 

получим  (индекс  p  внизу  означает  производную  по  этой  переменной),  что  искомые  функции  должны  иметь  вид  [Степанова,  2009  а;  2009  б]: 

$\begin{array}{c}{{\hat{\rho }}_1}^{\left(a\right)}(\omega ,q)={\rho }_1(\omega ,q,\lambda ),{{\hat{\rho }}_2}^{\left(a\right)}(\omega ,q,\lambda )=\text{}{\widehat{\stackrel{~}{}\rho }}_2(\omega ,q,\lambda ),\\ {\widehat{\stackrel{~}{}\rho }}_1(\omega ,q,\lambda )=-\frac{1}{2{\pi }^{}}\sum _{i=1}^N{\lambda }_i\underset{-\infty }{\overset{\infty }{\int }}\frac{{({\hat{Q}}_1^{\left(i\right)}{)}^{\text{'}}}_p^{}(\omega ,p-q)}{p}dp,\\ {\widehat{\stackrel{~}{}\rho }}_2(\omega ,q,\lambda )=-\frac{1}{2{\pi }^{}}\sum _{i=1}^N{\lambda }_i\underset{-\infty }{\overset{\infty }{\int }}\frac{{({\hat{Q}}_2^{\left(i\right)}{)}^{\text{'}}}_p^{}(\omega ,p-q)}{p}dp.\end{array}$ (9)

Таким  образом,  мы  приходим  к  следующей  системе  линейных  алгебраических  уравнений  (СЛАУ):

$A\lambda =f_{\delta },\lambda =({\lambda }_1,...,{\lambda }_N),f_{\delta }=(f_{1,\delta },...,f_{N,\delta }),$ (10) 

элементы  матрицы  которой  зависят  от ${\hat{Q}}_1^{\left(i\right)},{\hat{Q}}_2^{\left(i\right)}$ и ${\hat{Q}}_1^{\left(j\right)},{\hat{Q}}_2^{\left(j\right)}$.

3. АЛГОРИТМЫ ПОСТРОЕНИЯ АНАЛИТИЧЕСКИХ АППРОКСИМАЦИЙ АНОМАЛЬНЫХ ГРАВИТАЦИОННЫХ И МАГНИТНЫХ ПОЛЕЙ ИНТЕГРАЛОМ ФУРЬЕ (F-АППРОКСИМАЦИЯ)

В 50‑е–80‑е гг. в гравиметрии и магнитометрии благодаря исследованиям отечественных (К.В. Гладкий, Ф.М. Гольцман, В.Н. Луговенко, В.Р. Мелихов, С.А. Серкеров, Т.Н. Симоненко, О.А. Соловьев, В.Н. Страхов, Г.А. Трошков, С.В. Шалаев и др.) и зарубежных (B. Bhattacharryya, N. Barber, W. Dean, W. Hahn, P. Naidu, N. Neidell, F. Syberg и др.) ученых достаточно широко использовались методы обработки и интерпретации данных гравитационных и магнитных наблюдений, основанные на использовании аппарата спектрального анализа (анализа Фурье) элементов аномальных потенциальных полей.

При этом предполагалось, что элемент $V\left(x\right),x=(x_1,x_2,x_3)$ аномального поля непрерывно задан на всей бесконечной плоскости $x_3=0$ и что однозначно восстанавливается преобразование Фурье $F(u,v)$ элемента ${\left.V\left(x\right)\right|}_{x_3=0}:$

$F(u,v)=\frac{1}{2\pi }\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}{\left.V\left(x\right)\right|}_{x_3=0}\exp \left(i\right(u{x}_1+v{x}_2\left)\right)d{x}_{1}d{x}_2$  (11)

Спектральный анализ в гравиметрии и магнитометрии использовался для решения целого ряда интерпретационных задач. До самого последнего времени использование достижений математики шло путем простого заимствования, в подавляющем большинстве случаев путем использования готовых рецептов и аналитических конструкций, идей и методов. В результате теория и методы интерпретации, в том числе основанные на спектральном анализе геофизических данных, оказались неадекватными в реальной геофизической практике. Именно в существующей теории основными являются нереализуемые в вычислительной практике бесконечномерные задачи типа вычисления интегралов Фурье от аномальных полей.

В связи со сказанным естественно возникает вопрос – может быть, вообще методы, основанные на использовании спектральных представлений, т. е. на формуле (11) и обратной к ней:

$T\left\{V\left(x\right)\right\}=\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}K(u,v)F(u,v)\exp (-i(u{x}_1\text{}+v{x}_2\left)\right)dudv,$ (12)

где $T\left\{V\left(x\right)\right\}=W\left(x\right)$ есть некоторая линейная трансформанта функции $V\left(x\right)$, которой в спектральной области соответствует умножение спектра $F(u,v)$ на частотную характеристику $K(u,v),$ в реальных условиях несостоятельны? Ответ на этот вопрос отрицательный. А именно, если использовать тот общий подход, который был предложен В.Н. Страховым и назван им аппроксимационным, то ситуация радикальным образом изменяется. Создание численных методов нахождения спектров Фурье $F(u,v),$ основанных на общей теории метода линейных интегральных представлений, а также разработанных В.Н. Страховым теории и методах нахождения устойчивых приближенных решений линейных алгебраических уравнений больших размерностей, позволяет принципиально по-новому подойти к использованию метода анализа Фурье в задачах гравиметрии и магнитометрии [Страхов и др., 2009; Керимов, 2011]. Здесь важно подчеркнуть, что формулы (21) и (22) имеют абсолютно аналогичный вид и для четырехмерного пространства, частотная характеристика и спектр будут иметь еще одну компоненту:

$K\left(\omega \right)=K(u,v,w),F=F(u,v,w).$

Рассмотрим основные постановки задач на нахождение спектров Фурье элементов аномальных потенциальных полей по данным экспериментальных исследований этих полей. Для определенности ограничимся случаем гравитационного поля и задания значений одного элемента:

$\Delta g\left(x\right)=-\frac{\partial V_a\left(x\right)}{\partial x_3}$,  (13)

где $V_a\left(x\right),x=(x_1,x_2,x_3)$ суть потенциал аномального гравитационного поля, а ось $O{x}_3$ направлена вверх, в силу чего в (13) фигурирует знак минус.

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

Рассмотрим один из вариантов поставки F-апроксимации [Керимов, 2011], который состоит в том, что вводится спектральное представление функции $\frac{\partial V_a\left(x\right)}{\partial x_3}$, гармонической в полупространстве x3>–H, через спектр Фурье $F(u,v)$ потенциала V_a(x):

$\begin{array}{c}\frac{\partial V_a\left(x\right)}{\partial x_3}=\mathrm{Re}\left\{\frac{1}{2\pi }\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}K(u,v;x_3+H)F(u,v)\text{exp}(-i(u{x}_1+v{x}_2\left)\right)dudv\right\}=\\ =\frac{1}{2\pi }\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}K(u,v;x_3+H)\left(A\right(u,v\left)\text{cos}\right(u{x}_1+v{x}_2)+B(u,v\left)\text{sin}\right(u{x}_1+v{x}_2\left)\right)dudv\text{.}\end{array}$ (14)

В (14) положено:

$\begin{array}{c}K(u,v;x_3+H)=\\ =\sqrt{u^2+v^2}\exp \left(-(x_3+H)\sqrt{u^2+v^2}\right)\end{array}$ (15)

и

$F(u,v)=A(u,v)+iB(u,v).$ (16)

Основная вариационная постановка на нахождение функций A(u, v) и B(u, v) (действительной и мнимой частей комплексного спектра Фурье) и имеет следующий вид:

$\begin{array}{c}\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}{\left|F(u,v)\right|}^{2}dudv=\\ =\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}\left(A^2\right(u,v)+B^2(u,v\left)\right)dudv=\begin{array}{c}\\ \min \\ \begin{array}{c}A(u,v)\\ B(u,v)\end{array}\end{array}\end{array}$ (17)

при линейных условиях

$\begin{array}{c}{f}_{i,\delta }-\frac{1}{2\pi }\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}K(u,v;x_3^{\left(i\right)}+H)\times \\ \times \left[A\right(u,v\left)\text{\hspace{0.05em}\hspace{0.05em}cos}\right(u{x}_1^{\left(i\right)}+v{x}_2^{\left(i\right)})+\\ +B(u,v)\text{\hspace{0.05em}\hspace{0.05em}sin}(u{x}_1^{\left(i\right)}+v{x}_2^{\left(i\right)})]dudv=0\text{\hspace{0.17em}\hspace{0.17em},}\\ i=\text{1, 2, }...,N\text{,}\end{array}$ (18)

где положено

$f_{i,\delta }=V\left(x^{\left(i\right)}\right)+\delta V_i.$ (19)

Задача (17)–(19) решается методом множителей Лагранжа [Лаврентьев и др., 1950].

Значения параметров ${\lambda }_i$ (множителей Лагранжа) находятся из решения (СЛАУ) вида (10), в которой А есть (NxN) – матрица со свой­ством

$A=A^T\ge 0$ (20)

и элементами $a_{pq},$  $p=1,2,...,N,q=1,2,...,N:$

$a_{pq}\text{}=\text{}\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}\left[P_p\text{}(u,\text{}v)P_q(u,v)+\text{}{Q}_p(u,\text{}v)Q_q\text{}(u,v)\right]dudv.$  (21)

Можно показать [Страхов и др., 2009], что элементы матрицы (21) имеют вид:

 $a_{p,q}=\frac{6{(x_3^{\left(p\right)}+x_3^{\left(q\right)}+2H)}^3-9({(x_1^{\left(p\right)}-x_1^{\left(q\right)})}^2+{(x_2^{\left(p\right)}-x_2^{\left(q\right)})}^2)\cdot (x_3^{\left(p\right)}+x_3^{\left(q\right)}+2H)}{2\pi {\left[{(x_3^{\left(p\right)}+x_3^{\left(q\right)}+2H)}^2+{(x_1^{\left(p\right)}-x_1^{\left(q\right)})}^2+{(x_2^{\left(p\right)}-x_2^{\left(q\right)})}^2\right]}^{\frac{7}{2}}}$ .  (22)

Обратим внимание на тот факт, что условная экстремальная задача (20)–(22) не содержит априорной информации о свойствах погрешностей $\delta V_i$ в экспериментальных данных – в значениях $f_{i,\delta }$. Однако эта информация может быть учтена при нахождении устойчивых приближенных решений СЛАУ (10), но среднее значение помехи должно быть равно нулю.

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

$\begin{array}{c}\overline{f}(p,\phi )=\frac{1}{2\pi }\underset{-\infty }{\overset{+\infty }{\int }}F(u,v)e^{-i\left|\omega \right|p}d\left(\left|\stackrel{\rightharpoonup }{\omega }\right|\right),\stackrel{\rightharpoonup }{\omega }=(\cos \phi ,\sin \phi ),\\ F(u,v)=\underset{-\infty }{\overset{+\infty }{\int }}\overline{f}(p,\phi )e^{i\left|\omega \right|p}dp,\end{array}$ (23)

где через $\overline{f}(p,\phi )$ обозначено преобразование Радона от соответствующих переменных, а через  – преобразование Фурье. Если ввести обозначения и принять во внимание тот факт, что $\frac{ux+vy}{\sqrt{u^2+v^2}}=p,\phi =\arccos \frac{u}{\sqrt{u^2+v^2}}$, $\underset{-\infty }{\overset{+\infty }{\int }}{e}^{i\omega x}d\omega =2\pi \delta \left(x\right)$

(здесь $\delta \left(x\right)$ -дельта-функция Дирака), то можно показать, что:

 

 $\begin{array}{c}\underset{-\infty }{\overset{+\infty }{\int }}{e}^{-i(p-q){\rho }_1}\text{}{e}^{-iq{\rho }_2}dq\text{}=e^{-i{\rho }_{1}p}\underset{-\infty }{\overset{+\infty }{\int }}{e}^{i({\rho }_1-{\rho }_2)q}dq=\\ \text{}=\text{\hspace{0.17em}}2\pi \delta ({\rho }_1\text{}-{\rho }_2)e^{ip{\rho }_1},\\ 2\pi \underset{-\infty }{\overset{+\infty }{\int }}{e}^{ip{\rho }_1}\text{}{e}^{-{\rho }_1(z_i+H)}\text{}{e}^{i{\rho }_1(x_i\cos \phi +y_i\sin \phi )}d{\rho }_1\times \\ \times \underset{-\infty }{\overset{+\infty }{\int }}{e}^{-{\rho }_2(z_j+H)}{e}^{i{\rho }_2(x_j\cos \phi +y_j\sin \phi )}\delta ({\rho }_1-{\rho }_2)d{\rho }_2=\\ =\frac{2\pi }{{\left(\text{}{z}_i\text{}+2H\text{}+z_j\text{}\right)}^2\text{}+\text{}{\left(p\text{}-\left(\right(x_i\text{}+\text{}{x}_j)\cos \phi \text{}+(y_i\text{}+\text{}{y}_j\left)\sin \phi \right)\right)}^2}.\end{array}$ (24)

Если затем осуществить обратное преобразо­вание Радона, отнесенное к i-ой точке, то мы убедимся в том, что элементы матрицы системы линейных алгебраических уравнений (10) – это вторые производные по i-ой координате функции, преобразование Радона которой, деленное на $z_i+2H+z_j,$ есть как раз (24). Данное утверждение отражает очень важный факт: интегральные представления аномальных потенциальных полей (т. е. гармонических в некоторых областях пространст­ва истокообразно представимых функций) весьма тесно связаны друг с другом. Функция, описывающая рельеф поверхности, является непрерывной – поэтому ее можно рассматривать как предел некоторой гармонической в области, ограниченной этой поверхностью, функции. Если вспомнить выражения для элементов матрицы в методе S-аппроксимаций [Страхов и др., 2002а]:

$\begin{array}{c}{a}_{ij}=2\pi \left\{\frac{z_i+z_j}{{\rho }_{i,j}^3}+\frac{(z_i+z_j)(9{\rho }_{i,j}^2-6{(z_i+z_j)}^2)}{{\rho }_{i,j}^7}\right\},\\ {\rho }_{i,j}^2=(z_i+z_j{)}^2+(x_i-x_j{)}^2+(y_i-y_j{)}^2 1\le i,\\ j\le N,\end{array}$ (25)

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

Необходимо отметить, что не все постановки в методе F-аппроксимаций приводят к сформулированному выше утверждению (см., например, [Керимов, 2011]).

Фундаментальное решение уравнения Лапласа в  (равно как и производные этого решения по трем декартовым координатам) есть, как известно, свертка фундаментального решения в $R^3$ по четвертой координате:

$\begin{array}{c}\frac{z-x_3}{{\left({\left({(x-x_1)}^2+(y-x_2\right)}^2+{(z-x_3)}^2\right)}^{3/2}}=\\ =-\frac{1}{\pi }\frac{\partial }{\partial z}\left[\underset{-\infty }{\overset{+\infty }{\int }},\frac{dt}{{\left({(x\text{}-\text{}{x}_1)}^2\text{}+(y\text{}-\text{}{x}_2\right)}^2\text{}+{(z\text{}-\text{}{x}_3)}^2\text{}+{(t-x_4)}^2}\right].\end{array}$ (26)

Элементы матрицы системы типа (10) приобретают тогда вид (мы ограничились аппроксимацией только с помощью простого слоя):

${\tilde{a}}_{ij}^{}\text{}=\frac{4\pi (t_i+t_j+2H_t)}{{\left({(x_i\text{}-\text{}{x}_j)}^2\text{}+\text{}{(y_i\text{}-y_j)}^2\text{}+\text{}{(z_i\text{}-z_j)}^2\text{}+{(t_i\text{}+t_j\text{}+2H_t)}^2\right)}^2}.$  (27)

Поэтому любую функцию, заданную в некоторой области трехмерного пространства, можно рассматривать как след некоторой функции в четырехмерном пространстве при некотором значении четвертой координаты: $t=t_i,i=1,2,...,N.$ Такое рассмотрение позволяет нам находить значения гармонической в четырехмерном полупространстве функции, если задано распределение этой функции во всем трехмерном пространстве. Все утверждения, касающиеся связи различных интегральных преобразований в трехмерном декартовом пространстве, справедливы и для четырехмерного варианта благодаря инвариантности преобразований Радона и Фурье относительно параллельных переносов и вращений в декартовом пространстве любой размерности, большей или равной трем.

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

$\stackrel{⌢}{A}\left(x\right)\stackrel{⌢}{\lambda }=\stackrel{⌢}{f}\left(x\right),$ (28)

где матрица $\stackrel{⌢}{A}=\left\{{\stackrel{⌢}{a}}_{ij}\right\}$, $\stackrel{⌢}{\lambda }$ – решение системы линейных алгебраических уравнений в четырехмерном пространстве. Зависимость элементов матрицы от координат точки, в которой мы вычисляем значение поля, явная: она определяется по формуле (27).

Мы также применили подход, основанный на редукции свертки фундаментального решения в бесконечных пределах по четвертой координате к задаче с конечными пределами интегрирования, как это отмечал для двумерного преобразования Фурье В.Н. Страхов (см., напр., [Страхов и др., 2009]).

В этом случае коэффициенты матрицы (27) переходят в:

${\tilde{a}}_{ij}^{\left(2\right)}\text{}=\text{}\sum _{l=1}^L\frac{4\pi t_l\cdot \Delta t}{{\left({(x_i\text{}-\text{}{x}_j)}^2\text{}+\text{}{(y_i\text{}-\text{}{y}_j)}^2\text{\hspace{0.05em}}+\text{\hspace{0.05em}}{\left(t_l\right)}^2\text{}+{(z_i+z_j+2H_z)}^2\text{}\right)}^2}.$ (29).

Здесь $t_l,l=1,...,L,$ – некоторое множество точек на отрезке «вспомогательной оси» t. Шаг $∆t$ выбирается в зависимости от размера рассматриваемого полигона. При стремлении координаты t к бесконечности, формула (29) должна переходить в (25) (точнее, в ее первое слагаемое, соответствующее простому слою). Например, при размере участка, на котором задано поле, 60 км, желательно суммировать по $t_l$ от $-t_{\max }=-300$ до $t_{\max }$=300 км. Т. е. шаг по вспомогательной координате будет:

$\Delta t=\frac{2t_{\max }}{L}$.

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

Таким образом, мы приходим к идее метода интерпретации геофизических данных, основанного на попеременном переходе от одного пространства к другому: от «обычного» пространства трех измерений к объемлющему пространству R⁴ (или C²). Возникает вопрос как определить с приемлемой точностью значения поля в области, лежащей ниже источников? Сами источники в трехмерном пространстве относятся либо к области сингулярности поля (если речь идет о дискретных носителях), либо (в случае непрерывных распределений масс) мы должны рассматривать уравнение Пуассона.

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

1) сначала с помощью модифицированного метода S- и F-аппроксимаций находится эквивалентное распределение источников на некотором двумерном носителе (или на дискретном множестве носителей). Здесь важно подчеркнуть, что носитель необязательно должен представлять собой горизонтальную плоскость. Это может быть прямой двугранный угол (или двугранный угол иной радианной меры), часть сферической поверхности и т. п.;

2) затем мы находим распределение поля в области, внешней по отношению к носителю распределенных масс, но только в той части этой области, где достигается заданная интерпретатором относительная точность восстановления поля (например, не более 10 процентов);

3) по найденному (уже трехмерному!) распределению масс определяем распределение масс в четырехмерном пространстве в соответствии с формулой (49);

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

Шаги 1–4 можно повторять в итерационном цикле, задавая допустимые уровни погреш­ности.

4. КОМПЬЮТЕРНЫЕ ТЕХНОЛОГИИ НАХОЖДЕНИЯ ЛИНЕЙНЫХ АНАЛИТИЧЕСКИХ АППРОКСИМАЦИЙ ГАРМОНИЧЕСКИХ ФУНКЦИЙ (ЭЛЕМЕНТОВ ПОТЕНЦИАЛЬНЫХ ПОЛЕЙ) В ЛОКАЛЬНОМ ВАРИАНТЕ

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

4.1. Первый этап – формирование элементов матрицы А

На этом этапе формируется матрица А. Предварительно с использованием программ сортировки и выборки из исходных пунктов исключается определенное количество (N_контр) пунктов наблюдений. Конкретное количество и расположение исключенных пунктов для модельных и практических примеров будет указано при рассмотрении конкретных результатов ниже.

F-аппроксимация. Для F-аппроксимации расчетная формула для вычисления элементов матрицы А представляется в трехмерном пространстве следующим образом:

 $\begin{array}{c}{a}_{p,q}=\frac{3(2z_{p,q}^3-3{\rho }_{p,q}^2{z}_{p,q})}{2\pi {(z_{p,q}^2+{\rho }_{p,q}^2)}^{3.5}},z_{p,q}=x_3^{\left(p\right)}+x_3^{\left(q\right)}+\\ +2H,{\rho }_{p,q}=\sqrt{{(x_1^{\left(p\right)}-x_1^{\left(q\right)})}^2+{(x_2^{\left(p\right)}-x_2^{\left(q\right)})}^2}.\end{array}$ (30)

S-аппроксимация. Для S-аппроксимации расчетная формула для вычисления элементов матрицы А имеет вид (25).

По доказанному выше, абсолютно такие же формулы справедливы и для R-аппроксимации. В работах [Раевский и др., 2015а; 2015б] нами был предложен модифицированный метод S-аппроксимаций, который позволяет находить более устойчивые приближенные решения обратных задач геофизики и геоморфологии. Этот метод применялся при проведении расчетов, описанных в настоящей работе.

4.2. Решение СЛАУ

Решение СЛАУ является основной вычислительной проблемой при построении аппроксимационных конструкций на основе линейных интегральных представлений. В настоящей статье для решения СЛАУ используется регуляризованный итерационный трехслойный метод Чебышева с применением блочного метода контрастирования. Описание метода можно найти в работах [Степанова и др., 2016; Раевский, 2016].

4.3. Восстановление поля и нахождение его трансформант

F-аппроксимация. На данном этапе используя значения действительной A (u, v) и мнимой B (u, v) частей комплексного спектра Фурье, которые определяются с использованием формул (22) и значений компонент  вектора  множителей Лагранжа, вычисленных путем решения СЛАУ по описанным выше программам, вычисляются значения поля и его трансформант.

Рабочие формулы для вычисления восстановленного поля и нахождения его трансформант можно получить из формулы спектрального представления гравитационного поля в виде функции V(x), гармонической в полупространстве $x_3>-H$, через спектр Фурье:

$\begin{array}{c}V\left(x\right)=\frac{1}{2\pi }\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}\sqrt{u^2+v^2}{e}^{-z\sqrt{u^2+v^2}}\times \\ \times \left(A\right(u,v\left)\text{\hspace{0.05em}\hspace{0.05em}cos}\right(u{x}_1+v{x}_2)+B(u,v\left)\text{\hspace{0.05em}\hspace{0.05em}sin}\right(u{x}_1+v{x}_2\left)\right)dudv.\end{array}$ (31)

Подставив в данную формулу выражения для действительной A(u, v) и мнимой B(u, v) частей комплексного спектра Фурье, после очевидных преобразований получим соответствующие формулы.

Для вычисления первой вертикальной производной гравитационного потенциала используется следующая формула:

$V_z({\xi }_{1,l},{\xi }_{2,l},{\xi }_{3,l})=\frac{3}{2\pi }\sum _{k=1}^N{\lambda }_k\frac{(2z_{k,l}^2-3p_{k,l}^2)z_{k,l}}{{(z_{k,l}^2+p_{k,l}^2)}^{3.5}}$,  (32)

где:

$\begin{array}{c}{z}_{k,l}=x_{3,k}+{\xi }_{3,l}+2H,p_{k,l}=\\ =\sqrt{{(x_{1,k}-{\xi }_{1,l})}^2+{(x_{2,k}-{\xi }_{2,l})}^2},\end{array}$ (33) 

где: ${\xi }_l=({\xi }_{1,l},{\xi }_{2,l},{\xi }_{3,l})$ $(l=1,2,\dots ,M)$ - координаты пунктов, в которых вычисляются значения поля и его производных, $x_k=(x_{1,k},x_{2,k},x_{3,k})$ $(k=1,2,..,N)$ – координаты исходных пунктов наблюдений.

Покажем теперь, как комбинированный метод работает на модельных и практических примерах при решении обратных задач гравиметрии и геоморфологии [Тихонов и др., 1990; Ягола и др., 2014] (в частности, при построении аналитических аппроксимаций высотных отметок геоида).

4.4. Модельный пример № 2 по построению аналитический аппроксимации рельефа совместной S-, F-аппроксимацией

Функция $\Delta {\xi }_3\left({\xi }_1,{\xi }_2\right)$, описывающая рельеф земной поверхности, была условно представлена (т. е. рельеф так был смоделирован) как значения вертикальной производной потенциала большой группы (7000) точечных источников с различной избыточной плотностью, залегающих на расстоянии от 3.7 до 5.7 км ниже дневного рельефа (рис. 1), в точках горизонтальной плоскости $x_3=0$ км на равномерной сети точек $x^{\left(i\right)}=(x_1^{\left(i\right)},x_2^{\left(i\right)}),1\le i\le N.$ Общее количество точек измерений $N=9656,$ сеть нерегулярная. Из-за большого количества точечных источников смоделированный рельеф напоминает реальный в горной местности (рис. 2).

 

Рис. 1. Модельный пример № 2. Карта изолиний рельефа.

 

Рис. 2. Трехмерное изображение рельефа.

 

Представим функцию, описывающую рельеф, в виде суммы потенциалов простого и двойного слоев, распределенных на плоскости $x_3=-H$.

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

$\begin{array}{c}{a}_{ij}=\alpha a_{ij}^S+\left(1-\alpha \right)a_{ij}^F,\\ a_{ij}^S=2\pi \sum _{r=1}^r\left\{\frac{z_i+z_j-2H_r^S}{{\rho }_{i,j}^3}+\frac{(z_i+z_j-2H_r^S)(9{\rho }_{i,j}^2-6{(z_i+z_j-2H_r^S)}^2)}{{\rho }_{i,j}^7}\right\}\text{, }\\ {\rho }_{i,j}^2=(z_i+z_j+2H_r^S{)}^2+(x_i-x_j{)}^2+(y_i-y_j{)}^2 1\le i,a_{ij}^F=\sum _{k=1}^K\frac{3(2z_{i,j}^3-3{\rho }_{i,j}^2{z}_{i,j})}{2\pi {(z_{i,j}^2+{\rho }_{i,j}^2)}^{3.5}},\\ z_{i,j}=z_i+z_j-2H_k^F,{\rho }_{p,q}=\sqrt{{(x_i-x_j)}^2+{(y_i-y_j)}^2},\end{array}$ (34)

где: $0\le \alpha \le 1$ – весовой коэффициет; $a_{ij}^S$ – элементы матрицы для S-аппроксимации; $a_{ij}^F$ – элементы матрицы для F-аппроксимации; $M=\left(x_i,y_i,z_i\right),1\le i\le N$ – точки рельефа поверхности; $H_r^S,1\le r\le R$ – глубина залегания r-го носителя простого и двойного слоя в методе S-аппроксимаций; $H_k^F,1\le k\le K$ – характеристика полупространства представимости интеграла Фурье (то есть полагается, что поле создается совокупностью функций, гармонических при $z>H_k^F).$ Лучше полагать $R=K$ и $H_k^F=H_r^S,$ чтобы рельеф представлять в виде функций, каждая из которых гармоническая в определенной области. В дальнейшем будем считать, что $R=K$ и  $H_k^F=H_r^S=H_r.$

Для данного модельного примера проведены соответствующие вычисления при различных параметрах модели $(\alpha ,R,H_r).$ СЛАУ решалась при помощи регуляризованного итерационного метода Чебышева (CH). Мы провели также вычисления для данного примера по формуле (50) для L = 100. Результаты трехмерного и четырехмерного вариантов различаются приблизительно не более чем на 15% (в зависимости от разбиения вспомогательной оси t и сложности аппроксимируемого элемента поля или рельефа). В табл. 1 и табл. 2 содержатся результаты расчетов для различных значений L в четырехмерном варианте комбинированного метода S-, R- и F-аппроксимаций. Мы выполнили расчеты для тех же примеров (как модельных, так и практического), что и ранее для трехмерного варианта [Степанова и др., 2017], для того чтобы проиллюстрировать работу метода в четырехмерном пространстве. Таким образом, мы можем строить семейство методов аппроксимации данных детальных грави- и магниторазведки. Указанное семейство методов зависит от нескольких параметров: от коэффициента вхождения способа аппроксимации элемента аномального поля или функции, описывающей рельеф местности (S-, R- или F-аппроксимация), от размерности пространства, в котором решается обратная задача, от шагов разбиения по вспомогательным координатам.

 

Таблица 1. Результаты совместной S-F-аппроксимации для различных параметров модели. Модельный пример № 2

N N1 N2	$\sigma$min, м	$\sigma$max, м	Mодель	$\sigma$0, м	$\sigma$1, sконт1, м	$\sigma$2, sконт2, м	$∆$H, м	$\Sigma$/t, L = 150
9656 7832 8600	0.0322	0.4071	Mod_1	0.401	0.535 0.5746	0.4918 0.5218	4.0573	4.27339 · 10–4 8:17:09
			Mod_2	0.3218	0.499 0.5437	0.4536 0.4983	3.3829	3.4143 · 10–4 0:54:17
			Mod_3	0.2198	0.3782 0.4395	0.349 0.3899	2.1958	2.61289 · 10–4 0:48:55
			Mod_4	0.2196	0.3779 0.4392	0.3489 0.3897	2.1951	2.2378 · 10–4 0:48:57
			Mod_5	0.2195	0.3779 0.4391	0.3489 0.3896	2.1945	2.1561 · 10–4 0:47:49
			Mod_6	0.2196	0.3779 0.4391	0.3489 0.3896	2.1949	2.3568 · 10–4 0:48:58

Примечания: N – число исходных точек, N1 – число точек без контрольных точек I типа, N2 – число точек без контрольных точек II типа, ${\sigma }_1$ – среднеквадратичное отклонение, полученное с помощью аппроксимации по точкам рельефа без контрольных точек I типа; ${\sigma }_{конт1}$ – среднеквадратичное отклонение, полученное в контрольных точках I типа; ${\sigma }_2$ – среднеквадратичное отклонение, полученное с помощью аппроксимации по точкам рельефа без контрольных точек II типа; ${\sigma }_{конт2}$ – среднеквадратичное отклонение, полученное в контрольных точках II типа; $∆$H – максимальное расхождение между аппроксимированным и теоретически заданным значениями.

 

Таблица 2. Результаты совместной S-F-аппроксимации. Практический пример.Четырехмерный вариант

N N1 N2	$\sigma$min, м	$\sigma$max, м	Mодель	$\sigma$0, м	$\sigma$1, sконт1, м	$\sigma$2, sконт2, м	$∆$H, м	$\Sigma$/t, L = 200
10201 8038 9174	0.01	0.02	Mod_1	0.186	0.601 0.5794	0.5339 0.4972	1.473	5.03994 · 10–3 6:34:18
			Mod_2	0.198	0.612 0.5842	0.5437 0.5068	1.574	5.2531 · 10–3, 5:44:17

Примечания: N – число исходных точек, N1 – число точек без контрольных точек I типа, N2 – число точек без контрольных точек II типа; ${\sigma }_1$ – среднеквадратичное отклонение, полученное с помощью аппроксимации по точкам рельефа без контрольных точек I типа; ${\sigma }_{конт1}$ – среднеквадратичное отклонение, полученное в контрольных точках I типа; ${\sigma }_2$ – среднеквадратичное отклонение, полученное с помощью аппроксимации по точкам рельефа без контрольных точек II типа; ${\sigma }_{конт2}$ – среднеквадратичное отклонение, полученное в контрольных точках II типа; $∆$H – максимальное расхождение между аппроксимированным и теоретически заданным значениями.

 

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

Как видно из полученных результатов, наилучшие аппроксимации при выборе модели с несколькими плоскостями, расстояние между которыми относительно невелико (модель Mod_3). При этом результаты при различном значении весового коэффициента  мало отличаются друг от друга (см. результаты для Mod_3, Mod_4, Mod_5 и Mod_6). Вектор-решения, полученные для последних четырех моделей, почти не отличаются друг от друга – относительная норма разности менее 0.02, а коэффициент корреляции более 0.99, а значит решения можно считать почти тождественными. Поэтому в дальнейшем мы будем использовать модель с весовым коэффициентом   0.5.

4.5. Построение аналитической аппроксимации участка горного рельефа

Данные по рельефу для апробации алгоритмов по аппроксимации рельефа предоставлены заказчиком. Исходные данные представляют собой цифровые модели рельефа (ЦМР) в системе координат WGS84 размером 30 x 30 км (или меньше, в зависимости от конкретной ЦМР) с шагом 5 м. ЦМР подразделяются на три типа, в зависимости от степени разреженности рельефа: горы, равнины и холмы. Для апробации алгоритма выбран горный участок местности размером 19 x 30 км. Количество точек по оси OX $N_x=3841$, а по оси OY – $N_y=6041$. Общее количество точек в ЦМР $N=N_x\cdot N_y=23203481$. ЦМР представляет из себя детальную горную цепь, максимальный перепад высот составляет более 3400 м (рис. 3). Благодаря столь крупному масштабу (1:500), на трехмерном изображении исходной ЦМР можно разглядеть даже отдельные тропинки в горах.

Исходная ЦМР слишком большая для построения аппроксимации на персональном компьютере, для апробации алгоритма из данной ЦМР выбрана область 500 x 500 м (от исходной ЦМР от 0 до 500 м по обеим осям) с шагом 5 м, общее число точек $N_s=10201$. Были построены соответствующие аппроксимации рельефа поверхности.

4.6. Пример построения модифицированных S-аппроксимаций данных гравиметрической съемки

В данном примере исходные данные аномального гравитационного поля Земли (рис. 4) в области от 0 до 55$°$ северной широты и от 120 до 180$°$ восточной долготы характеризуются сложным строением поля, амплитуда поля изменяется от –15 до 80 мГал. Сетка регулярная, шаг сетки 0.5$°$, общее число пунктов измерений N=13431. Рассматриваемая область соответствует северной части Тихого океана, охватывающая наиболее глубокие океанические желоба (марианский, алеутский, курило-камчатский, японский, филиппинский).

При построении S-аппроксимаций использовался способ трехступенчатого контроля, который был предложен В.Н. Страховым. В качестве априорной информации были выбраны следующие значения для оценки квадрата нормы ошибки ${\delta }_{\min }^2=300,{\delta }_{\max }^2=2000.$ В рамках метода S-аппроксимации аномальное гравитационное поле аппроксимировалось суммой простого и двойного слоев, распределенных на сфере радиуса $R_0=6365$ км. Результаты модифицированных S-аппроксимаций показаны в табл. 3.

Как видно, использование блочного метода контрастирования сократило время вычислений в 1,5 раза. Сравним полученные решения. Коэффициент корреляции ${||{\tilde{x}}_1-{\tilde{x}}_2||}_E/\left({||{\tilde{x}}_1||}_E+\left({||{\tilde{x}}_2||}_E\right)\right)=0.985$ и отношение норм ${||{\tilde{x}}_1||}_E/{||{\tilde{x}}_2||}_E=1.20$ (${\tilde{x}}_1$ – решение, полученное методом CH, ${\tilde{x}}_2$ – решение методом CH + блочный метод). Заметим, что значение среднеквадратического отклонения  при решении только методом CH значительно выше, чем при решении методом CH + блочный метод. В силу этого для дальнейшего анализа решение, полученное при использовании блочного метода контрастирования, будет приоритетнее, чем решение, полученное решением СЛАУ только итерационным методом.

Для выявления региональных особенностей, построено аналитическое продолжение поля на R6367 км и карта первой вертикальной производной.

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

 

Рис. 3. Карта изолиний рельефа выбранного участка горного рельефа.

Рис. 4. Карта наблюденного аномального гравитационного поля северной части Тихого океана.

 

Таблица 3. Результаты S-аппроксимации

Метод	$\sigma$min, мГал	$\sigma$max, мГал	$\sigma$1, sконтр1, мГал, N1 = 10 456	$\sigma$2, sконтр2, мГал, N2 = 11 872	$\sigma$0, мГал, N = 13 431	Δ, t
CH	0.1495	0.3859	0.3771 0.4189	0.3673 0.3914	0.3556	1.0856 · 10–2, 31:12
CH+блочный метод	0.1495	0.3859	0.3061 0.3324	0.2978 0.3122	0.2910	8.8826 · 10–3, 20:47

Примечания: $\sigma$₀ – среднеквадратичное отклонение, полученное с помощью аппроксимации по исходным точкам, N – общее количество точек измерений; ${\sigma }_1$ – среднеквадратичное отклонение, полученное с помощью аппроксимации по точкам измерений без контрольных точек первого типа; $\sigma$_контр1 – среднеквадратичное отклонение, полученное в контрольных точках первого типа; N1 – общее количество точек измерений без учета контрольных точек первого типа; $\sigma$₂ – среднеквадратичное отклонение, полученное с помощью аппроксимации по точкам измерений без контрольных точек второго типа; $\sigma$_контр2 – среднеквадратичное отклонение, полученное в контрольных точках второго типа; N2 – общее количество точек измерений без учета контрольных точек второго типа. Расхождение между $\sigma$₀ и $\sigma$_контр1, $\sigma$_контр2 составляет менее 1.20 при решении методом CH и менее 1.15 при решении методом CH + блочный метод, что свидетельствует о высокой качественности аппроксимации в минимальных по модулю точках исследуемого поля.

 

5. ЗАКЛЮЧЕНИЕ

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

 

Рис. 5. Линейные трансформанты поля: (а) – первая вертикальная производная uz; (б) – аналитическое продолжение поля на R=6367 км. Белым выделены наиболее контрастные области, являющиеся зонами субдукций.

 

Из полученных результатов по построению аналитической аппроксимации рельефа (разделы 4.4., 4.5.) видно, что совместная S-F-аппроксимация позволяет построить аналитическую аппроксимацию рельефа горной поверхности с достаточно высокой степенью точности (среднеквадратическое отклонение ${\sigma }_0\approx 20$ см) (см. табл. 1, табл. 2). Построение аналитической аппроксимации горного рельефа – одна из самых трудных задач, поэтому общее вычислительное время здесь достаточно велико (см. табл. 2).

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

Авторы благодарят РФФИ за финансовую поддержку данной работы (грант 17-51-53002‑ГФЕН-а).

 

About the authors

I. E. Stepanova

Institute of the Earth Physics of the Russian Academy of Sciences

Author for correspondence.
Email: tet@ifz.ru
Russian Federation, Bolshaya Gruzinskaya str., 10-1, Moscow 123242, Russia

I. A. Kerimov

Institute of the Earth Physics of the Russian Academy of Sciences

Email: tet@ifz.ru
Russian Federation, Bolshaya Gruzinskaya str., 10-1, Moscow 123242, Russia

A. G. Yagola

Lomonosov Moscow State University

Email: tet@ifz.ru
Russian Federation, 1, Leninskie gory, Moscow, 119991

References

Supplementary files