Физика Земли, 2021, № 4, стр. 148-157

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

И. Э. Степанова^1, 2, *, И. А. Керимов¹, А. А. Спесивцев^1, 2, 3, В. А. Тимофеева^1, 2, 4, П. С. Михайлов^1, 2

¹ Институт физики Земли им. О.Ю. Шмидта РАН
г. Москва, Россия

² Научно-технологический университет “Сириус”
г. Сочи, Россия

³ ФГБУ “Центр геодезии, картографии и ИПД”
г. Москва, Россия

⁴ Институт теории прогноза землетрясений и математической геофизики РАН
г. Москва, Россия

^* E-mail: tet@ifz.ru

Поступила в редакцию 17.01.2021
После доработки 19.02.2021
Принята к публикации 03.03.2021

DOI: 10.31857/S0002333721040098

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

Аннотация

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

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

ВВЕДЕНИЕ

При изучении гравитационного поля Земли и решении большого числа прикладных задач полезной оказывается информация о топографии в силу известной корреляции аномалий силы тяжести и функции, описывающей рельеф поверхности или грид глубин [Конешов, Проценко 1984].

Целью настоящей статьи является выполнение спектрального анализа профильных топографических данных, а также модельных данных о рельефе на том же участке, полученных с помощью глобальной модели гравитационного поля EIGEM-6C4 и построенных на основе спутниковых данных цифровых моделей рельефа (ЦМР). При совпадении спектральных характеристик реальных топографических данных (измеренных с помощью аппаратуры, установленной на различных транспортных средствах: автомобилях, самолетах, морских судах и т.п.) и значений высотных отметок рельефа поверхности Земли или грида глубин Мирового океана можно упростить решение задач геофизики и геодезии интерпретационного характера. Профильные данные обладают, если можно так выразиться, весьма “ограниченными возможностями” с точки зрения применения их при исследованиях глубинного строения Земли и поиске месторождений полезных ископаемых. Для интерпретации данных гравимагниторазведки гораздо более предпочтительными являются так называемые площадные данные, но проведение подобного рода съемки физических полей нашей планеты – дело весьма затратное. Поэтому профильные данные о поле или рельефе можно попытаться дополнить значениями, полученными с помощью спутниковых систем наблюдения и (или) вычисленными в соответствии с формулами, принятыми в той или иной модели гравитационного поля Земли.

Различные варианты модифицированных S- и F-аппроксимаций (см. [Страхов, Керимов 1999; Степанова и др., 2019; 2020б]) позволяют построить линейные аналитические аппроксимации элементов аномальных потенциальных полей и функции, описывающей рельеф местности.

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

Нами рассматривались задачи с большим количеством точек (десятки тысяч). Возникающие системы линейных алгебраических уравнений (СЛАУ) решались с помощью усовершенствованного метода блочного контрастирования [Степанова, 2019]

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

МЕТОДИКА ИССЛЕДОВАНИЯ

Алгоритмы построения аналитических аппроксимаций с помощью интеграла Фурье (F-аппроксимация)

Аппарат спектрального анализа (анализа Фурье) элементов аномальных потенциальных полей довольно широко применялся в 50–80-х годах прошлого века благодаря исследованиям как отечественных, так и зарубежных ученых (см., например, [Осколков, 1997; Болотин, Вязьмин, 2018; Зорич, 1984; Осипов, 2019]).

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

(1)

$F(u,v) = \frac{1}{{2\pi }}\int\limits_{ - \infty }^{ + \infty } {\int\limits_{ - \infty }^{ + \infty } {{{{\left. {V(x)} \right|}}_{{{{x}_{3}} = 0}}}} } \exp (i(u{{x}_{1}} + v{{x}_{2}}))d{{x}_{1}}d{{x}_{2}}.$

Обратное преобразование имеет вид:

(2)

$\begin{gathered} T\{ V(x)\} = \frac{1}{{2\pi }} \times \\ \times \,\,\int\limits_{ - \infty }^{ + \infty } {\int\limits_{ - \infty }^{ + \infty } {K(u,v)F(u,v)\exp ( - i(u{{x}_{1}} + v{{x}_{2}}))dudv} } , \\ \end{gathered} $

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

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

(3)

$\Delta g(x) = - \frac{{\partial {{V}_{a}}(x)}}{{\partial {{x}_{3}}}}.$

Здесь ${{V}_{a}}\left( x \right),x = \left( {{{x}_{1}},{{x}_{2}},{{x}_{3}}} \right)$ суть потенциал аномального гравитационного поля, а ось $0{{x}_{3}}$ направлена вверх (в “воздух”), в силу чего в (3) фигурирует знак минус. Будем считать, что рельеф (высотные отметки или грид глубин) представляет собой предельные значения некоторой гармонической в нижнем (соответственно – верхнем) полупространстве функции, которую будем обозначать так же, как элемент гравитационного поля. Подчеркнем тот факт, что задача определения спектра элемента потенциального поля или функции, описывающей рельеф, рассматривается в трехмерном пространстве. Таким образом, речь идет о вычислении двумерного спектра полезного сигнала.

Постановка задачи состоит в том, что вводится спектральное представление функции $\frac{{\partial {{V}_{a}}(x)}}{{\partial {{x}_{3}}}}$, гармонической в полупространстве x₃ > –H или ${{x}_{3}} < H$ через спектр Фурье $F(u,v)$ потенциала V_a(x):

(4)

$\begin{gathered} \frac{{\partial {{V}_{a}}(x)}}{{\partial {{x}_{3}}}} = \operatorname{Re} \left\{ {\frac{1}{{2\pi }}\int\limits_{ - \infty }^{ + \infty } {\int\limits_{ - \infty }^{ + \infty } {K(u,v;{{x}_{3}} + H)F(u,v){\text{exp}}( - i(u{{x}_{1}} + v{{x}_{2}}))dudv} } } \right\} = \\ = \frac{1}{{2\pi }}\int\limits_{ - \infty }^{ + \infty } {\int\limits_{ - \infty }^{ + \infty } {K(u,v;{{x}_{3}} + H)(A(u,v){\text{cos}}(u{{x}_{1}} + v{{x}_{2}}) + B(u,v){\text{sin}}(u{{x}_{1}} + v{{x}_{2}}))dudv} } {\text{.}} \\ \end{gathered} $

Здесь положено:

(5)

$K(u,v;{{x}_{3}} + H) = \exp ( - ({{x}_{3}} + H)\sqrt {{{u}^{2}} + {{v}^{2}}} )$

(6)

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

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

(7)

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

(8)

$\begin{gathered} {{f}_{{i,\delta }}} - \frac{1}{{2\pi }}\int\limits_{ - \infty }^{ + \infty } {\int\limits_{ - \infty }^{ + \infty } {K(u,v;x_{3}^{{(i)}} + H)} } \times \\ \times \,\,[A(u,v){\text{cos}}(ux_{1}^{{(i)}} + vx_{2}^{{(i)}}) + \\ + \,\,B(u,v){\text{sin}}(ux_{1}^{{(i)}} + vx_{2}^{{(i)}})]{\kern 1pt} dudv = 0{\text{,}} \\ i = 0{\text{,1,2,}}...{\text{,}}N{\text{,}} \\ \end{gathered} $

при этом компоненты сигнала имеют вид:

(9)

${{f}_{{i,\delta }}} = V({{x}^{{(i)}}}) + \delta {{V}_{i}}.$

Задача (7)–(8) решается методом множителей Лагранжа. Двумерные спектры полезного сигнала можно представить в виде:

(10)

$A(u,v) = \sum\limits_{i = 1}^N {{{\lambda }_{i}}} {{P}_{i}}(u,v){\text{,}}\,\,\,\,B(u,v) = \sum\limits_{i = 1}^N {{{\lambda }_{i}}{{Q}_{i}}(u,v){\text{,}}} $

где положено:

(11)

${{P}_{i}}(u,v) = \frac{1}{{2\pi }}K(u,v,x_{3}^{{(i)}} + H){\kern 1pt} {\text{cos}}{\kern 1pt} (ux_{1}^{{(i)}} + vx_{2}^{{(i)}}),$

(12)

$\begin{gathered} {{Q}_{i}}(u,v) = \frac{1}{{2\pi }}K(u,v,x_{3}^{{(i)}} + H){\kern 1pt} {\text{sin}}{\kern 1pt} (ux_{1}^{{(i)}} + vx_{2}^{{(i)}}){\text{,}} \\ i = 1{\text{,2,}}...{\text{,}}N{\text{.}} \\ \end{gathered} $

С учетом (5) имеем:

(13)

$\begin{gathered} A(u,v) = \frac{1}{{2\pi }}\sum\limits_{i = 1}^N {{{\lambda }_{i}}{{e}^{{ - (x_{3}^{{(i)}} + H)\sqrt {{{u}^{2}} + {{v}^{2}}} }}}\cos (ux_{1}^{{(i)}} + vx_{2}^{{(i)}})} = \\ = \sum\limits_{i = 1}^N {{{\lambda }_{i}}{{P}_{i}}} (u,v), \\ \end{gathered} $

(14)

$\begin{gathered} B(u,v) = \frac{1}{{2\pi }}\sum\limits_{i = 1}^N {{{\lambda }_{i}}{{e}^{{ - (x_{3}^{{(i)}} + H)\sqrt {{{u}^{2}} + {{v}^{2}}} }}}\sin (ux_{1}^{{(i)}} + vx_{2}^{{(i)}})} = \\ = \sum\limits_{i = 1}^N {{{\lambda }_{i}}{{Q}_{i}}} (u,v), \\ \end{gathered} $

где:

(15)

${{P}_{i}}(u,v) = \frac{1}{{2\pi }}{{e}^{{ - (x_{3}^{{(i)}} + H)\sqrt {{{u}^{2}} + {{v}^{2}}} }}}\cos (ux_{1}^{{(i)}} + vx_{2}^{{(i)}}),$

(16)

${{Q}_{i}}(u,v) = \frac{1}{{2\pi }}{{e}^{{ - (x_{3}^{{(i)}} + H)\sqrt {{{u}^{2}} + {{v}^{2}}} }}}\sin (ux_{1}^{{(i)}} + vx_{2}^{{(i)}}).$

Как показано в работе [Страхов, Керимов, 1999], значения параметров ${{\lambda }_{i}}$ (множителей Лагранжа) находятся из решения системы линейных алгебраических уравнений (СЛАУ):

(17)

${\text{A}}\lambda = {{f}_{\delta }},$

в которой А есть (N × N) – матрица со свойством:

(18)

${\text{A}} = {{{\text{A}}}^{{\text{T}}}} \geqslant 0$

и элементами ${{a}_{{pq}}},$ $1 \leqslant p{\text{, }}q \geqslant N{\text{:}}$

(19)

${{a}_{{pq}}} = \int\limits_{ - \infty }^{ + \infty } {\int\limits_{ - \infty }^{ + \infty } {\left[ {{{P}_{p}}(u,v){{P}_{q}}(u,v) + {{Q}_{p}}(u,v){{Q}_{q}}(u,v)} \right]} dudv} ,$

а ${{f}_{\delta }}$ есть N – вектор с компонентами ${{f}_{{i,\delta }}}$, определенными по (9); $\lambda $ – есть N – вектор с компонентами ${{\lambda }_{i}}$.

В работе [Страхов, Керимов 1999] были получены аналитические выражения для элементов матрицы (19):

(20)

${{a}_{{pq}}} = \frac{1}{{4\pi {{{\left( {x_{3}^{{\left( p \right)}} + x_{3}^{{\left( q \right)}} + 2H} \right)}}^{3}}{{{\left[ {{{{\left( {x_{3}^{{\left( p \right)}} + x_{3}^{{\left( q \right)}} + 2H} \right)}}^{2}} + {{{\left( {x_{1}^{{\left( p \right)}} - x_{1}^{{\left( q \right)}}} \right)}}^{2}} + {{{\left( {x_{2}^{{\left( p \right)}} - x_{2}^{{\left( q \right)}}} \right)}}^{2}}} \right]}}^{{\frac{1}{2}}}}}}.$

Решая систему (17), мы находим множители Лагранжа ${{\lambda }_{i}},\,\,i = 1,2, \ldots ,N$ и, тем самым, спектры искомых элементов аномальных полей или рельефа $K(u,v,{{x}_{3}} + H)F(u,v),$ ${{\left| {K(u,v,{{x}_{3}} + H)F(u,v)} \right|}^{2}}$. Последнее выражение представляет собой зависимость квадрата амплитуды спектра от частот u и v – энергетический спектр на высоте H (или глубине Н, если рассматривается грид глубин). Именно нахождение квадрата амплитуды спектра различных типов сигналов и позволяет впоследствии графически изобразить функцию плотности вероятности.

При применении структурно-параметрического подхода функции $\rho _{r}^{{}}(\xi )$, возникающие в рамках вариационного подхода, претерпевают некоторые изменения. Именно, соотношение:

(21)

$\tilde {\rho }_{r}^{{}}(\xi {\text{;}}\lambda ) = {{p}_{r}}(\xi )\sum\limits_{i = 1}^N {{{\lambda }_{i}}Q_{r}^{{(i)}}(\xi )} ,\,\,\,\,r = 1,2,...,R$

просто обобщается следующим образом:

(22)

Смысл конструкции (22) состоит в том, что параметры, от которых зависят аналитические выражения искомых функций, делаются зависящими от индекса r.

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

(23)

в которой $\mathop \lambda \limits^ \circ $ есть блочный вектор вида:

(24)

а матрица ${\AA}$ – это блочная матрица следующего вида:

(25)

${\AA} = \begin{array}{*{20}{c}} \vline & {\,\,{{{\AA}}^{{(1)}}}\,\,}&\vline & {\,\,{{{\AA}}^{{(2)}}}\,\,}&\vline & {\begin{array}{*{20}{c}} {_{{_{{_{{_{{_{{}}}}}}}}}}}&{^{{^{{^{{^{{}}}}}}}}} \end{array}}&\vline & {\,\,{{{\AA}}^{{(R)}}}\,\,}\vline & \end{array},$

в которой блоки ${{{\AA}}^{{(r)}}}$ имеют элементы:

(26)

и обладают свойством:

(27)

${{{\AA}}^{{(r)}}} = {{{\AA}}^{{(r),T}}} \geqslant 0,\,\,\,\,r = 1,2, \ldots ,R.$

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

Алгоритмы нахождения аналитических аппроксимаций профильных топографических данных в локальном варианте

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

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

На этом этапе формируется матрица А. Предварительно с использованием программ сортировки и выборки из исходных пунктов исключается определенное количество (N_контр) пунктов наблюдений. Профильные топографические данные характеризуются тем, что на отрезке (профиле), имеющем достаточно большую для локального варианта длину, имеется значительное число отметок рельефа (в нашем случае число данных наблюдений превышало 40 000). Мы разбили исходный трек (профиль) на несколько участков таким образом, чтобы в каждом блоке при применении усовершенствованного метода блочного контрастирования [Степанова и др., 2020б] было не более 10 000 строк. Такой подход позволяет эффективно использовать оперативную память одного персонального компьютера (процессора) или, в случае проведения параллельных вычислений, рационально распределять задачи между различными процессорами.

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

(28)

${{a}_{{pq}}} = \frac{1}{{4\pi {{z}^{3}}{{{\left( {z_{{pq}}^{2} + \rho _{{pq}}^{2}} \right)}}^{{\frac{1}{2}}}}}},$

где

$\begin{gathered} {{z}_{{pq}}} = x_{3}^{{\left( p \right)}} + x_{3}^{{\left( q \right)}} + 2H, \\ {{\rho }_{{pq}}} = \sqrt {{{{\left( {x_{1}^{{\left( p \right)}} - x_{1}^{{\left( q \right)}}} \right)}}^{2}} + {{{\left( {x_{2}^{{\left( p \right)}} - x_{2}^{{\left( q \right)}}} \right)}}^{2}}} . \\ \end{gathered} $,

Второй этап – нахождение устойчивого приближенного решения СЛАУ

Разработка алгоритмов получения устойчивых приближенных решений плохо обусловленных СЛАУ является ключевой проблемой при построении аналитических аппроксимаций в рамках метода линейных интегральных представлений. В.Н. Страховым была создана новая теория регуляризации СЛАУ, основные положения которой изложены в большой серии работ [Страхов и др., 2009; Страхов, Страхов, 1999а; 1999б; и др.]. Эти методы решения СЛАУ реализованы в пакете прикладных программ на языке С для СЛАУ с симметрической положительно полуопределенной матрицей и приближенно заданной правой частью (авторы – В.Н. Страхов, А.В. Страхов).

Для решения больших и сверхбольших систем линейных алгебраических уравнений нами был усовершенствован метод блочного контрастирования [Степанова и др., 2020а], основанный на предложенном В.Н. Страховым методе МБКоорС – метод блочного координатного спуска [Страхов, Страхов, 1999а].

Этот метод может быть использован во всех трех случаях: а) переопределенных систем (N > M); б) нормально определенных систем (N = M = n); в) недоопределенных систем (N < M).

Суть МБКоорС состоит в том, что матрица A, соответственно вектор x в исходной системе (17) представляются в блочной форме:

(29)

$A = \begin{array}{*{20}{c}} \vline & {\,\,{{{\text{A}}}_{1}}\,\,}&\vline & {\,\,{{{\text{A}}}_{2}}\,\,}&\vline & {\begin{array}{*{20}{c}} {\,\,\,\,\,\,...}&{{{{_{{_{{_{{_{{}}}}}}}}}}^{{^{{}}}}}} \end{array}}&\vline & {\,\,{{{\text{A}}}_{Q}}\,\,}\vline & \end{array},$

(30)

$x = \begin{array}{*{20}{c}} \vline & {{{x}^{{(1)}}}}\vline & \\ \hline \vline & {{{x}^{{(2)}}}}\vline & \\ \hline \vline & {}\vline & \\ \hline \vline & {{{x}^{{(Q)}}}}\vline & \end{array},$

где в (29)–(30) числа блоков Q удовлетворяют условию

(31)

$Q \ll M,$

а размерности блоков ${{{\text{A}}}_{q}}$ матрицы A и блоков ${{x}^{{(q)}}}$ вектора x, q = 1, 2, …, Q суть:

(32)

$(N \times {{M}_{q}}),\,\,{{M}_{q}}$

соответственно. В общем случае все ${{M}_{q}}$ могут быть различны, но важнейшими являются два частных случая:

(33)

$1)\,\,{{M}_{q}} = c = {\text{const}},\,\,\,\,q = 1,2, \ldots ,Q$

(34)

$2)\,\,{{M}_{q}} = c = {\text{const}},\,\,\,\,q = 1,2, \ldots ,Q - 1,\,\,\,\,{{M}_{Q}} < c;$

при этом всегда должны иметь место неравенства

(35)

$N > {{M}_{q}},\,\,\,\,q = 1,2, \ldots ,Q.$

В обоих случаях число c должно быть выбрано таким образом, чтобы блоки ${{{\text{A}}}_{q}}$ и векторы ${{f}_{\delta }}$ целиком размещались в оперативной памяти компьютера.

В усовершенствованном методе блочного контрастирования благодаря учету априорной информации мы можем получить решение, адекватное реальной геофизической практике. Фактически, мы решаем задачу нахождения устойчивого приближенного решения СЛАУ (17) с помощью трехслойного усовершенствованного чебышевского итерационного метода для каждого из блоков q = 1,2, …, Q либо с помощью регуляризации метода Холецкого. При этом размерность каждого из блоков меняется в соответствии с выбранным нами законом (линейным или, при большой размерности исходного массива данных о рельефе, нелинейным – как правило, степенным).

Третий этап – восстановление рельефа

На данном этапе, используя значения действительной A(u, v) и мнимой B(u, v) частей комплексного спектра Фурье, которые определяются с использованием формул (13)–(14) и значений компонент ${{\lambda }_{i}}$ вектора $\lambda $ множителей Лагранжа, вычисленных описанным выше путем решения СЛАУ, вычисляются значения рельефа.

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

(36)

$\begin{gathered} V(x) = \frac{1}{{2\pi }}\int\limits_{ - \infty }^{ + \infty } {\int\limits_{ - \infty }^{ + \infty } {{{e}^{{ - z\sqrt {{{u}^{2}} + {{v}^{2}}} }}}(A(u,v)\,{\text{cos}}\,(u{{x}_{1}} + v{{x}_{2}})} + } \\ + \,\,B(u,v)\,{\text{sin}}\,(u{{x}_{1}} + v{{x}_{2}}))dudv. \\ \end{gathered} $

Подставив в данную формулу выражения для действительной A(u, v) и мнимой B(u, v) частей комплексного спектра Фурье, после очевидных преобразований, получим соответствующие формулы. В формуле (36) речь идет о трехмерных полях и двумерных спектрах. Наши численные эксперименты показали, что предложенная методика модифицированных F-аппроксимаций на основе структурно-параметрического подхода работает и при анализе данных профильной съемки. Одна из частот в формуле (36), скажем, v, полагается равной нулю (поскольку у-компонента измеренного сигнала также равна нулю). Таким образом, мы получаем одномерный спектр функции, описывающей высотные отметки на профиле. Необходимо отметить, что отсутствие площадной съемки не привело к принципиальной невозможности спектрального анализа рельефных отметок с помощью метода F-аппроксимаций.

Спектральный анализ профильных топографических данных, полученных с помощью ГНСС, и сравнение спектров реального и модельного рельефов

Для определения спектральных характеристик и выявления общих составляющих в сигналах “разной” природы был выбран набор данных, составленный на основе измерений, полученных во время переезда автомобиля-лаборатории ИФЗ РАН по маршруту Москва–Беломорск (24–25.06.2019 г.).

В качестве основного измерительного оборудования использовался комплект высокоточной ГНСС-аппаратуры. Подробное описание комплекса приборов и их характеристики приводятся в работе [Передерин и др., 2018]. Для обеспечения лучших условий наблюдений ГНСС-антенна устанавливалась на крыше автомобиля-лаборатории. Во время движения по маршруту велась непрерывная запись измерений с дискретностью 0.1 с (10 Гц), что соответствует пространственному разрешению 2.5 м при скорости движения автомобиля 90 км/ч.

На основе полученных ГНСС-измерений были определены геодезические координаты (широта, долгота и высота) точек маршрута с использованием метода абсолютных точных координатных определений Precise Pont Positioning (PPP).

В качестве альтернативного источника информации о высотах точек земной поверхности использовались данные, полученные с использованием ЦМР ArcticDEM [Showstack, 2017]. Cпециально построенная на основе данных спутниковой миссии Sentinel-1 модель также может применяться при решении геофизических задач интерпретационного характера: авторами настоящей работы был поставлен ряд математических экспериментов по построению модифицированных S-аппроксимаций топографических данных, полученных с помощью указанной модели. Относительная точность аналитических аппроксимаций превысила требуемую на практике и составила ${{10}^{{ - 9}}}.$ Государственно-частный проект ArcticDEM Национального агентства геопространственной разведки (National Geospatial-Intelligence Agency, США) и Национального научного фонда (National Science Foundation, США) был реализован в 2016 г. для автоматического создания высококачественной ЦМР с высоким разрешением территории суши к северу от 60° северной широты на основе оптических стереоизображений посредством высокопроизводительных вычислений и программного обеспечения для фотограмметрии с открытым исходным кодом. Данные модели ArcticDEM предоставляются в виде сетки геодезических высот с пространственным разрешением от 2 до 1000 м. Значения высот вдоль точек маршрута определялись на основе билинейной интерполяции с использованием модуля grdtrack из пакета программ Generic Mapping Tools. Также для областей пролегания трека экспериментальных измерений ИФЗ РАН 2019 г. по маршруту Москва–Беломорск были построены ЦМР с помощью радиолокационной интерферометрии на основе использования разности фаз спутниковых снимков. Для работы использовались данные спутниковой миссии Sentinel-1 Европейского космического агентства. Отбор снимков для построения ЦМР производился с учетом сроков проведения измерений автомобилем-лабораторией (даты отобранных снимков – 14.06.2019 г. и 26.06.2019 г.), а также оптимальных значений временной (12 дней) и пространственной (76.03 м) базовых линий. Полученная ЦМР имеет геометрическое разрешение 5 × 20 м.

Данные об аномалиях силы тяжести определялись с использованием вычислительного сервиса Международного центра глобальных моделей Земли (ICGEM) на основе глобальной комбинированной модели геопотенциала EICGEM-6С4 представленной в виде сферических гармоник до 2190 степени. Значения аномалий силы тяжести определялась на сетке с шагом в 0.1 град. Значения вдоль точек маршрута определялись на основе билинейной интерполяции.

Таким образом, перечень анализируемых данных состоял из:

1) высот, определенных по данным ГНСС, в метрах;

2) высот, определенных по ЦМР ArcticDEM, в метрах;

3) высот, определенных по ЦМР, построенной на основе данных спутниковой миссии Sentinel-1, в метрах;

4) аномалий силы тяжести, определенных по модели ГПЗ EIGEM-6C4, в мГал.

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

$\begin{gathered} 1)\,\,{{M}_{q}}(t) = \alpha t + {{M}_{0}},{\kern 1pt} \,\,\alpha > 0,\,\,\,{{M}_{q}}(0) = {{M}_{0}}; \\ {{M}_{q}}(1) = {{M}_{1}}. \\ 2)\,\,{{M}_{q}}(t) = \beta \sqrt {t + {{t}_{0}}} ,\,\,\beta > 0,\,\,{{M}_{q}}(0) = {{M}_{0}};\,\,{{M}_{q}}(1) = {{M}_{1}}. \\ \end{gathered} $

Минимальный размер блока (соответствующий индексу “0”) выбирался равным 1000, максимальный (соответствующий индексу “1”) – 10 000.

F-аппроксимации топографических данных были выполнены со среднеквадратическим отклонением, равным 1.5–3 см, при этом максимальное различие между высотными отметками наблюдалось для профиля 1-го типа (полученного с помощью установленной на автомобиле системы ГНСС) и составило 5.6 см. Следует отметить, что программное обеспечение, которое было разработано нами для решения поставленной в данной работе задачи, позволяет аппроксимировать высотные отметки рельефа с большей точностью (до 0.001 см), однако на практике такой уровень погрешности задавать не следует: нужно учитывать реальные условия измерений.

Рис. 1.

Маршрут переезда Москва–Беломорск.

РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЯ

В результате были оценены амплитуды периодических составляющих (синфазных и квадратурных) в рядах высот, определенных на основе ГНСС-данных, для диапазона частот от 0 до 9 Гц (рис. 2–рис. 4).

Рис. 2.

График квадратов амплитуд косинусоидальной (синфазной) составляющей сигнала высот, определенных из данных ГНСС.

Рис. 3.

График квадратов амплитуд синусоидальной (квадратурной) составляющей сигнала высот, определенных из данных ГНСС.

Рис. 4.

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

При этом в спектре синфазной составляющей наибольшую амплитуду обнаружили низкочастотные гармоники, соответствующие примерно 0.1 и 0.8 Гц. Что при длине анализируемого участка маршрута в 500 км соответствует длинам волн в 5000 и 625 км.

В спектре квадратурной составляющей наибольшую амплитуду обнаружили низкочастотные гармоники, соответствующие примерно 0.2 и 0.6 Гц. Что при длине анализируемого участка маршрута в 500 км соответствует длинам волн в 2500 и 830 км.

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

Рис. 5.

Плотность вероятности спектра данных, определенных по ГНСС-измерениям.

Также был проведен расчет топографии по данным модели ArcticDEM (рис. 6), а также гравитационного поля (по EIGEM-6C4), заданного в тех же точках, в которых измерялись высотные отметки рельефа. Значимых различий в соответствующих спектрах не наблюдается.

Рис. 6.

Плотность вероятности спектра данных, определенных по модели ArcticDEM.

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

В статье приводится эффективный алгоритм спектрального анализа топографических профильных данных полученных с помощью ГНСС и ЦМР, основанных на разных типах данных (оптических и радиолокационных).

Описанная методика построения аналитических аппроксимаций позволяет сравнивать спектральные характеристики рельефа полученного на основе ГНСС-измерений с соответствующими модельными значениями. Причем существенных различий между использованием значений модели ArcticDEM и модели, построенной на основе данных спутниковой миссии Sentinel-1, не наблюдается. Что позволяет сделать вывод, что при подобных задачах на усмотрение исследователей может быть выбрана как модель ArcticDEM, так и специально построенная для конкретного участка рельефа и в соответствии с датами проведения измерений ГНСС ЦМР на основе данных спутниковых миссий, например имеющейся в открытом доступе миссии Европейского космического агентства Sentinel-1.

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

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

Болотин Ю.В., Вязьмин В.С. Спектральный анализ точности векторной аэрогравиметрии // Фундаментальная и прикладная математика. 2018. Т. 22. Вып. 2. С. 33–57.
Зорич В. А. Математический анализ. М.: Физматлит. 1984. 544 с.
Керимов И.А., Степанова И.Э., Раевский Д.Н. Комбинированные аппроксимационные методы решения задач гравиразведки и магниторазведки // Геология и геофизика Юга России. 2018. № 3. С. 37–50.
Конешов В.Н., Проценко С.В. Корреляционные связи геофизических полей. М.: ИФЗ АН СССР. 1984. 100 с.
Конешов В.Н., Степанова И.Э. Аналитические аппроксимации грида глубин арктического бассейна и их спектральный анализ // Физика Земли. 2008. № 5. С. 34–41.
Осипов О.В. Спектральный анализ дискретных сигналов с высоким частотным разрешением // Вычислительные методы и программирование. 2019. Т. 20. Вып. 3. С. 270–282.
Осколков К.И. Рельефная аппроксимация, анализ Фурье–Чебышева и оптимальные квадратурные формулы. Теория приближений. Гармонический анализ / Никольский С.М. Труды МИАН. 1997. Т. 219. С. 269–285.
Передерин Ф.В., Алешин И.М., Иванов С.Д., Михайлов П.С., Погорелов, В.В., Холодков К.И. Портативный комплекс регистрации сигналов ГНСС: полевые испытания и перспективы использования // Наука и технологические разработки. 2018. № 4. Т. 97. С. 28–40.
Степанова И.Э., Керимов И.А., Раевский Д.Н., Щепетилов А.В. Комбинированный метод F-, S-и R-аппроксимаций при решении задач геофизики и геоморфологии // Физика Земли. 2018а. № 1. С. 96–113.
Степанова И.Э., Керимов И.А., Раевский Д.Н., Щепетилов А.В. О совершенствовании методов обработки больших объемов данных в геофизике и геоморфологии на основе модифицированных S- и F-аппроксимаций // Физика Земли. 2020а. № 3. С. 82–97.
Степанова И.Э., Керимов И.А., Ягола А.Г. Аппроксимационный подход в различных модификациях метода линейных интегральных представлений // Физика Земли. 2019. № 2. С. 31–47.
Степанова И.Э., Раевский Д.Н., Щепетилов А.В. Эффективная технология построения цифровых моделей рельефа и аналитических аппроксимаций потенциальных полей Земли // Наука и технологические разработки. 2018б. Т. 97. № 3. С. 5–14.
Степанова И.Э., Щепетилов А.В., Погорелов В.В., Михайлов П.С. Структурно-параметрический подход при построении цифровых моделей рельефа и гравитационного поля Земли с использованием аналитических S-аппроксимаций // Геофизические процессы и биосфера. 2020б. Т. 19. № 2. С. 107–116.
Страхов В.Н., Керимов И.А. Аппроксимационная реализация спектрального анализа в гравиметрии и магнитометрии. Основные проблемы теории интерпретации гравитационных и магнитных полей / Страхов В.Н. М.: ОИФЗ РАН. 1999. С. 183–206.
Страхов В.Н., Керимов И.А., Степанова И.Э. Разработка теории и компьютерной технологии построения линейных аналитических аппроксимаций гравитационных и магнитных полей. М.: ИФЗ РАН. 2009. 254 с.
Страхов В.Н., Страхов А.В. О решении систем линейных алгебраических уравнений с приближенно заданной правой частью, возникающих при решении задач гравиметрии и магнитометрии // Изв. секции наук о Земле РАЕН. 1999а. Вып. 3, ноябрь. С. 20–49.
Страхов В.Н., Страхов А.В. Центральная вычислительная проблема в рамках третьей парадигмы в теории и практике интерпретации потенциальных полей. 1. Основы теории и численных методов. О некоторых вопросах теории интерпретации потенциальных полей. М.: ОИФЗ РАН, 1999б. С. 128–191.
Showstack R. Map provides high-resolution look at nearly entire Arctic region // Eos. 2017. V. 98. https://doi.org/10.1029/2017EO085451
Stepanova I.E., Kerimov I.A., Shchepetilov A.V. Improving the methods for processing large data in geophysics and geomorphologybased on the modified S- and F-approximations // Izv., Phys. Solid Earth. 2020. V. 56. № 3. P. 364–378.
Stepanova I.E., Raevskiy D.N., Shchepetilov A.V. Separation of potential fields generated by different sources // Izv., Phys. Solid Earth. 2020. V. 56. № 3. P. 379–391.

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

Инструменты

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

Физика Земли

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