Моделирование аномалий силы тяжести эквивалентными источниками в полярных областях Земли

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

В настоящее время 180 статических моделей геопотенциала Земли представлены на сайте Немецкого научно-исследовательского цен- тра наук о Земле в городе Потсдам (ICGEM), который является одним из шести центров Международной гравитационной полевой службы Международной ассоциации геодезии. Все эти модели доступны в виде нормированных коэффициентов сферических гармоник (Balmino et al., 2012) в стандартном формате, который был принят Европейским космическим агентством (ESA – European Space Agency) в качестве официального формата данных международного космического проекта GOCE (Global Ocean Circulation Experiment).

Одним из инструментов работы с глобальными моделями ГПЗ является выделение различных компонент поля в редукции Буге Д^ б (трансформация) и анализ их мор-

Работа лицензирована в соответствии с CC BY 4.0. Чтобы просмотреть копию

этой лицензии, посетите

фологии, направленный на выявление устойчивых взаимосвязей между аномалиями силы тяжести и особенностями геологического строения площади исследований (Петрищев-ский, 2023; Терентьев и др., 2024; Чадаев и др., 2023). Целью трансформации является преобразование цифровых моделей поля △ д_Б, направленное на «фокусировку» информации об аномалиях, связанных с отдельными геологическими телами, их группами или определенными геологическими границами (Гравиразведка: Справочник геофизика, 1990). Общепринятой классификации методов трансформации гравитационного и магнитного полей нет. Среди большого числа существующих линейных трансформаций часто выделяют две основные группы, различающиеся по своим спектральным характеристикам: 1) группу «региональных трансформаций», предназначенных для построения регионального фона и сглаживания исходного поля; 2) группу «локальных трансформаций», предназначенных для выделения локальных аномалий различных порядков, обусловленных сравнительно небольшими по размерам объектами, залегающими на малых глубинах (Долгаль, 2022).

Есть два класса вычислительных методов, использующихся для определения региональных и локальных компонент геопотенциальных полей: 1) класс, в котором задачи трансформации решаются на основе традиционных методов вычислительной математики (Фурье-преобразования и цифровой фильтрации сигналов); 2) класс, где те же задачи решаются на основе аппроксимации (приближения) полей потенциальными (истокообразными) функциями, представляющими собой эффекты эквивалентных источников (Гравиразведка: Справочник геофизика, 1990). Далее будет рассматриваться метод трансформации, который относится к классу 2, предложенному в 60-х гг. ХХ в. и успешно применяемому целым рядом исследователей (Долгаль, Пугин и др., 2022). Его суть заключается в приближении наблюденного поля △gs теоретическим полем Ддт, представляющим собой аномальный эффект системы элементарных тел, массы определяются путем решения системы линейных алгебраических уравнений (СЛАУ) с приближенно заданной правой частью:

Gm = AgB,           (1)

где G = {gij} — квадратная матрица значений гравитационных эффектов для точечного источника с единичной массой (т = 1); m -вектор неизвестных значений аномальных масс; △gs – вектор значений аномалий силы тяжести (в данном случает – 1-й радиальной производной гравитационного потенциала VR). Под гравитационным эффектом gij в данном случае подразумевается радиальная производная гравитационного потенциала точечного источника VR, определенная в сферической системе координат £<р,Л,г: dV/dR=VR(R0,pM = f (Ro-rcoswyr3

где f = 6,67х10-11 м3-кг1-с-2 - гравитационная постоянная; R0, ^0, Л0 - координаты точки измерений; r,^,A  - координаты источ ника; r0 = VRq + r2 — 2R0rcosw, ® - угол при центре O земного шара между точкой измерений и источником:    созш = cosy0cosp + sin^0sin^ cos(A0 — Л).

Полярные области Земли – это территории, расположенные на крайнем севере и крайнем юге нашей планеты и занимающие в совокупности примерно 1/12 ее поверхности. В пределах этих территорий возникают специфические сложности, связанные с решением СЛАУ (1) для GRID-моделей ГПЗ, характеризующихся регулярным шагом по широте B и долготе L . Особенности моделирования аномалий силы тяжести △ д Б системой точечных масс будут рассмотрены далее.

Особенности моделирования гравитационного поля в полярных областях

Истинная фигура Земли является геоидом, который при картографических построениях аппроксимируется эллипсоидом вращения. Геодезические координаты (широта B , долгота L ) относятся к общеземному эллипсоиду, размеры и форма которого определяются значениями большой полуоси и сжатия (для России – ПЗ-90.11). Экспериментально установлено, что степень отличия

«эллипсоидальной» и «сферической» моделей Земли при выполнении трансформаций гравитационного поля не превышает 0,05 % (Долгаль, Костицын и др., 2022). Таким образом, при региональных исследованиях можно ограничиться представлениями о ша-роообразной форме Земли. Для точек задания поля долгота L остается неизменной, но происходит замена геодезической широты B на геоцентрическую широту Ф. Для перехода от геодезической широты к геоцентрической целесообразно использовать формулу Каврайского: Ф = В — 8'39" sin2B (Белкин и др., 1988). В свою очередь, при решении СЛАУ (1) вместо геоцентрической широты используется коширота ^ = 90° — Ф.

Рис. 1. Вертикальные сечения эллипсоидальной и сферической моделей Земли (плоскость Гринвичского меридиана): геодезическая широта B (а) и геоцентрическая широта Ф (б): красный цвет – вектор силы тяжести

Обозначим шаг GRID-модели гравитационного поля в градусной мере по широте - А 1 , по долготе - А 2 и введем функцию r ( А ) для преобразования градусной меры в линейную. Известно, что хорошее качество аппроксимации аномалий силы тяжести и высокую скорость сходимости итерационного процесса решения СЛАУ (1) обеспечивает соблюдение условия 1 < R * /r(A) <1,5, где R* - удаление точечной массы от точки задания поля по радиусу Земли вниз; г( Л ) - расстояние между точками измерений: r(A) = r(A 1 ) = r(A 2 ) (Аронов, 1976). Для глобальных моделей ГПЗ Земли в условиях высоких широт при △ 1 = А 2 возникает неравенство r(A 1 ) ^ r(A 2 ), поэтому данное условие не выполняется (рис. 2). Это приводит к существенному увеличению чисел обусловленности cond G матрицы коэффициентов СЛАУ при R* = const и создает дополнительные вычислительные сложности.

Обусловленность матрицы коэффициентов СЛАУ является важнейшей характеристикой, определяющей сложность процесса ее решения и точность полученных результатов (Бахвалов и др., 2000). Число обусловленности ( Н -число Тодда) матрицы G :

cond G = vG = ||G-1|| • ||G|| (2)

определяет влияние неточностей 5А^В в задании исходных данных на окончательный результат:

||6m||/||m|| =vc HSAgBH/HAgBH. (3)

Напомним, что минимально возможная величина v_G = 1, СЛАУ с высокими значениями V g >> 1 называют плохо обусловленными. В пределе относительная погрешность решения может в V g раз превысить погрешность задания исходных данных, которая обычно составляет не менее первых единиц мГал (Михайлов и др., 2022), т.е. около 1-2 % .

Рис. 2. Изменение длины 1 ° дуги параллели от экватора к полюсу Земли

Рассмотрим изменения размеров площади модели поля и характеристик матрицы G при перемещении от экватора к Северному полюсу по меридиану. Размер модели в градусной мере составляет 20 ° ^х20 ° , шаг перемещения - 10 ° . Сеть точек задания поля -0,5 ° ^х0,5 ° , глубина эквивалентных источников (точечных масс) R* = г(Д 1 ) = 55 км.

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

2 0,5 Фробениуса ПСП / = pS /to/ ) J .

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

Таблица 1. Результаты вычислительного эксперимента: влияние удаления модели поля с фиксированным размером 20 ° ^х20 ° от экватора Земли

Местоположение и площадь модели			Характеристики матрицы коэффициентов
Южная граница	Северная граница	Площадь, кв. км	Ранг r(G)	Норма ПСП /	Число vG
10 ° ю.ш.	10 ° с.ш.	4948158	1681	0.497	21
0 °	20 ° с.ш.	4870506	1681	0.500	22
10 ° с.ш.	30 ° с.ш.	4640684	1681	0.509	26
20 ° с.ш.	40 ° с.ш.	4267428	1681	0.526	35
30 ° с.ш.	50 ° с.ш.	3764571	1681	0.556	64
40 ° с.ш.	60 ° с.ш.	3149814	1681	0.606	194
50 ° с.ш.	70 ° с.ш.	2443561	1681	0.689	2617
60 ° с.ш.	80 ° с.ш.	1667859	1502	0.849	2431481

Как видно, площадь модели при удалении от экватора закономерно сокращается почти в 3 раза, т.к. ее линейные размеры по долготе L уменьшаются. За счет сближения точек задания поля для полярной области Земли отмечается резкое увеличение чисел обусловленности матрицы коэффициентов (vG > 2,4 х 106), что влечет за собой вычислительные сложности при решении СЛАУ (1). При наиболее близком к полюсу местоположении модели число линейно независимых векторов-строк в системе (1) сокращается на ~10 %, ее ранг становится меньше числа неизвестных, что влечет за собой бесконечное множество возможных решений.

Моделирование гравитационного поля Сибирской платформы

При выполнении исследований авторами статьи использовалась одна из моделей ГПЗ – EIGEN-GRGS.RL04.MEAN-FIELD, полученная на основе данных спутниковых миссий GRACE и SLR в 2019 г.: гравитационные аномалии А^б в полной редукции Буге, определенные для плотности материков 2,67 г/см3 и плотности океанов 1,05 г/см3. Для территории размером ~5,15 млн км2, ограниченной 52°-72° с.ш., 84°-132° в.д., охватывающей всю Сибирскую платформу и структуры ее обрамления, были получены значения поля с шагом - 0,4°. Размерность GRID-модели – 51×121 (N = 6171 точка). Источником информации о рельефе земной поверхности являлась модель ETOPO1 (Amante, Eakins, 2009). Диапазон изменения аномалий силы тяжести ~387 мГал, высотные отметки H лежат в пределах от –18 до 2420 м (рис. 3).

Рис. 3. Карты изоаномал силы тяжести VR в полной редукции Буге (а) и изогипс рельефа H земной поверхности (б)

Глубина эквивалентных источников была принята равной шагу сети точек поля по меридиану R* = г(А1) = 55 км. Истокообразная аппроксимация и решение СЛАУ выполнялись тремя разными способами. Первый («традиционный») способ: размещение точечных масс по каждой точкой задания поля и решение СЛАУ из 6171 уравнения с 6171 неизвестным модифицированным методом Гаусса (решатель solve библиотеки Numpy языка Python (Хилл, 2021)). Второй способ: создание равномерной сети точечных масс с шагом г ^ 55 км и таким же значением R*, с последующим решением СЛАУ с 4681 неизвестным методом наименьших квадратов (решатель lstsq библиотеки Numpy языка Python (Хилл, 2021)). Третий способ: использование регуляризованного метода сингулярного разложения (процедура svd библиотеки Numpy языка Python (Хилл, 2021)) при рабо- те с аппроксимационной конструкцией, представленной в описании способа 1.

Сингулярное разложение (Singular Values Decomposition, SVD) показывает геометрическую структуру матрицы и позволяет наглядно представить имеющиеся данные (Васин, 2020), в частности оценить погрешности при работе с матрицами неполного ранга или близкими к вырождению. С помощью SVD можно определить тип и меру обусловленности матриц, вычислить обратную матрицу. Усеченное сингулярное разложение матриц применяется для решения плохо обусловленных СЛАУ (Андрушев-ский, 2008). Матрицу коэффициентов G можно представить в виде:

G = USVT, (4) где U, V - ортогональные матрицы (UT = U-1, VT = V-1); S - диагональная матрица с коэффициентами 8i (Si > <5i+1), которые называются сингулярными числами матрицы G. С использованием (4) решение СЛАУ (1) определяется следующим образом:

Число обусловленности матрицы G (2) можно найти с использованием минимального и максимального значений ее собственных чисел X (Бахвалов и др., 2000):

m = VS-1U-1△gБ.         (5)

^VG    ^ max /^ min .

Достаточно вычислять только r векторов-столбцов U^–1 и r векторов-строк V , соответствующих числу r наибольших сингулярных чисел S_r > S_m i _n, где S_m i _n - выбранное пороговое значение. Согласно теореме Эккарта-Янга, таким образом получается наилучшее приближение G⁺ к матрице G^-1 на множестве всех матриц ранга, не превосходящего r.

С другой стороны, собственные числа матрицы G^TG являются квадратами сингулярных чисел, полученных при разложении (4), что позволяет использовать соотношение A i = S i при оценке V g и выборе r . В данном случае было принято r = 5500 (рис. 4).

номер строки матрицы S

Рис. 4. Сингулярные числа матрицы коэффициентов G при моделировании аномалий силы тяжести Сибирской платформы

Информация по решению СЛАУ описанными выше способами представлена в табл. 2. Наименьшее по невязке решение имеет способ 1, однако он характеризуется наибольшим по норме результатом m . Результаты решения систем способами 2 и 3 характеризуются среднеквадратическим отклонением (СКО), близким к предполагаемой точности использованной модели ГПЗ, значения ||m|| y для них на порядок ниже.

Результаты пересчета гравитационного поля в верхнее полупространство на высоту

100 км (рис. 5а, в) для аналитических моделей 2 и 3 практически совпадают. Плохая обусловленность СЛАУ и отсутствие регуляризации приводят к заметным искажениям результата пересчета поля на высоту для модели 1 (рис. 5б). Для изображения фрагмента территории (рис. 5б, в), ограниченного координатами 66,8 ° -72 ° с.ш., 84 ° -98 в.д., использовалась равнопромежуточная цилиндрическая проекция, отражающая непосредственно цифровые данные, без сглаживающего влияния интерполяции.

Таблица 2. Информация по решению СЛАУ при моделировании аномалий силы тяжести Сибирской платформы

Способ	V g	Невязка, мГал				IlmH /	Время счета, с
		Минимум	Максимум	Среднее	СКО
1	1155694	–0,04	0,03	0	0,001	1,84×107	269
2	220900	–38,36	41,72	0	4,15	0,78×106	223
3	250000	–28,32	29,17	0	2,88	1,17×106	432

Рис. 5. Результаты пересчета гравитационного поля в верхнее полупространство на высоту 100 км: для всей территории способом 2 (а) и для ее фрагмента (красный контур) способами 1 (б) и 3 (в)

Заключение

Применение метода эквивалентного источника к глобальным GRID-моделям ГПЗ в пределах полярных областей Земли осложнено плохой обусловленностью СЛАУ вида (1), решение которых требуется для определения численных значений аппроксимирующих масс. Это является следствием резкого уменьшения расстояния r(A 2 ) между точками задания поля по долготе по мере увеличения широты B . Простым и достаточно эффективным способом борьбы с ним является ранее предложенное в работе (Долгаль, Новикова и др., 2022) уменьшение глубин R ^* источников поля с ростом B . Однако его использование приводит к появлению высокочастотных компонент в спектре аналитической модели поля вблизи полюсов Земли.

Неустойчивость решений СЛАУ при условии размещения источников под каждой точкой поля и R* = const продемонстрирована в данной статье на примере Сибирской платформы. Однако можно построить близ- кую к регулярной в плане сеть размещения эквивалентных источников, шаг которой будет соответствовать их глубине R* = r(L1) = г(Д2). Затем для решения нормальной СЛАУ GTGm = GT△gБ применить метод наименьших квадратов. При данном способе аппроксимации размерность модели источников всегда будет меньше размерности цифровой модели ГПЗ.

Также впервые апробирован подход, базирующийся на сингулярном разложении матрицы коэффициентов G, отвечающей «традиционной» аппроксимационной конструкции. Метод SVD, который называют «томографом высокого разрешения», входит в десятку наиболее важных алгоритмических достижений ХХ в. (Андрушевский, 2008). Он позволяет определить ранг матрицы G, меру ее обусловленности, вычислить обобщенную обратную матрицу G+, найти оптимальное решение несовместной системы уравнений с матрицами неполного ранга, а также осуществить приближенное решение систем уравнений с плохо обусловленной матрицей пол- ного ранга. По сравнению с процедурой linalg.lstsq сингулярное разложение linalg.svd требует повышенных в 1,5–1,8 раза затрат машинного времени. Но наглядная возможность формирования матрицы G+ с требующимся числом обусловленности vG, по мнению авторов, полностью компенсирует эти затраты.

Практическое применение представленных алгоритмов – изучение глубинного строения, структурно-тектоническое районирование и прогнозирование полезных ископаемых в Арктике и Антарктике на основе информации, представленной в спутниковых (альтиметри-ческих) глобальных моделях ГПЗ. Следует добавить, что на практике для истокообразной аппроксимации необходимо использовать минимум две телескопированные цифровые модели гравитационного поля с разным шагом сети при различных уровнях размещения точечных масс (Долгаль, Рыжов, 2024).

Исследование выполнено при финансовой поддержке Министерства науки и высшего образования РФ в рамках государственного задания (рег. номер НИОКТР: 124020500054-3).