Применение матричной декомпозиции для построения аналитических моделей аномалий силы тяжести

Долгаль Александр Сергеевич

doi:doi:

Главная / Журналы / Вестник АГГИ / Применение матричной декомпозиции для построения аналитических моделей аномалий силы тяжести

Отправить рукопись Скачать PDF
Текст

Цитировать

Цитирований:

Применение матричной декомпозиции для построения аналитических моделей аномалий силы тяжести

Журнал: ВЕСТНИК АГГИ

Рубрики: ИНЖЕНЕРНАЯ И РУДНАЯ ГЕОФИЗИКА

Долгаль Александр Сергеевич ¹

Информация об авторах и публикации

Авторы:

1. Горный институт Уральского отделения РАН (EGPГлавный научный сотрудник)
сотрудник с 01.01.2002 по настоящее время

Уральский государственный горный университет

Свердловская область, Россия

Тип:

Статья

Страницы:

с 1 по 1

Статус:

Классификаторы:

ВАК 1.6 Науки о Земле и окружающей среде
ВАК 1.2.1 Искусственный интеллект и машинное обучение
ВАК 1.2.2 Математическое моделирование, численные методы и комплексы программ
УДК 55 Геология. Геологические и геофизические науки
УДК 550 Вспомогательные геологические науки
УДК 550.3 Геофизика

Язык материала:

русский

Ключевые слова:

гравитационное поле, аномалии Буге, аналитическая модель поля, точечные массы, система линейных уравнений, разложение Холецкого, сингулярное разложение, QR-разложение, регуляризация, Арктика, Антарктика.

Финансирование:

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

Аннотация и ключевые слова

Аннотация:
Рассмотрено применение матричной декомпозиции для построения аналитических моделей региональных аномалий силы тяжести с помощью точечных масс на сферообразной Земле. Размещение точечных масс осуществляется на двух уровнях глубин. Разложение Холецкого, SVD и QR-разложения эффективно используются при решении систем линейных алгебраических уравнений в процессе истокообразной аппроксимации. Представлен пример моделирования гравитационного поля в редукции Буге для территории Сибирской платформы размером около 5 млн. кв. км с координатами 52°–72° с. ш., 84°–132° в. д.

Ключевые слова:
гравитационное поле, аномалии Буге, аналитическая модель поля, точечные массы, система линейных уравнений, разложение Холецкого, сингулярное разложение, QR-разложение, регуляризация, Арктика, Антарктика.

Текст

Текст (PDF): Читать Скачать

Введение

Под аналитической моделью аномалий силы тяжести $∆g$ в данном случае подразумевается множество параметров сеточной модели эквивалентных источников (точечных масс), порождающих модельное гравитационное поле $∆g_T$ , весьма близкое наблюденному полю $‖∆g-∆g_T ‖≤ε$ , где величина $ε$ определяется предполагаемой нормой помех $δ$ , $δ=‖ε‖$ . В это множество входят геометрические и физические параметры, т. е. априорно заданные координаты размещения источников и значения их масс соответственно. Будем рассматривать моделирование аномалий силы тяжести в редукции Буге на фрагменте сферообразной Земли. По классификации В. Н. Страхова, речь пойдет о региональном варианте линейных аналитических аппроксимаций, который применяется для сравнительно больших территорий размером 10⁵–10⁶ км² и более [18]. В этом случае в качестве модели Земли используется сфера Каврайского, для расчетов выполняется переход от геодезических координат (широты B и долготы L) к сферическим координатам $∑R$ , а аномалии силы тяжести $∆g$ отождествляются с 1-й радиальной производной гравитационного потенциала $V_R$ [7].

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

Gm = V_R, (1)

где G — квадратная матрица значений гравитационных эффектов для точечного источника с единичной массой $(m=1)$ , V_R — вектор наблюденных аномалий силы тяжести.

Если при решении практических задач порядок СЛАУ (1) составляет N ≌ (1 – 5) ´ 10⁴ и менее, то для ее решения можно отказаться от использования приближенных итерационных методов [1] и перейти к алгоритмам, основанным на матричной декомпозиции [13]. Этот вариант достаточно часто возникает при работе с глобальными моделями гравитационного поля Земли [12, 14], разрешающая способность которых ограничивает число независимых значений аномалий силы тяжести в пределах изучаемых территорий. Далее будет рассмотрено применение разложения Холецкого, QR-разложения и сингулярного разложения (SVD — singular values decomposition) в процессе истокообразной аппроксимации региональных гравитационных аномалий $V_R$ в редукции Буге. Суть данного подхода заключается в представлении матрицы G в виде произведения двух и более матриц, обладающих определенными свойствами (ортогональностью, симметричностью, диагональностью), облегчающими решение задачи (1).

Методы разложения матриц

Методы матричной декомпозиции, по мнению, высказанному в журнале SIAM Rewiev в 2000 г., входят в список из 10 алгоритмов, оказавших наибольшее влияние на науку и индустрию в ХХ веке. Фундаментальные работы по этой проблеме были опубликованы учеными из США Полом Дуайером (Paul Dwyer) и Альстоном Скоттом Хаусхолдером (Alston Scott Householder) в 1951 и 1954 гг. Далее рассмотрим 3 способа факторизации матриц: разложение Холецкого, QR-разложение и сингулярное разложение.

Разложение Холецкого позволяет представить положительно определенную симметричную матрицу коэффициентов нормальной системы уравнений в виде произведения левой и правой треугольных матриц G^TG = LL^T, где L — нижняя треугольная матрица со строго положительными элементами главной диагонали. QR-разложение — это представление невырожденной матрицы G в виде G = QR, где Q — ортогональная матрица (Q^T = Q^-1), a R — верхняя треугольная матрица [13]. QR-разложение выполняется с сохранением нормы Фробениуса исходной матрицы. Это является его преимуществом по сравнению с классическим подходом в методе наименьших квадратов, где требуется решение нормальной СЛАУ вида G^TGm = G^TV_R. Число обусловленности матрицы коэффициентов такой системы достаточно велико: cond G^TG ≈ cond² G .

Сингулярное разложение (Singular Values Decomposition, SVD) позволяет представить матрицу G размером $k×k$ в виде, представленном на рисунке 1.

G = USV^T, (2)

где U, V — ортогональные матрицы (U^T = U^-1, V^T = V^-1), $(U^T=U^{-1},V^T=V^{-1})$ S $=\{δ_{ii} \}$ — диагональная матрица с элементами $δ_i (δ_i≥δ_{i+1})$ , которые называются сингулярными числами матрицы [13].

Первый в нашей стране опыт применения сингулярного разложения для факторного анализа петрофизических данных и решения линейной обратной задачи гравиразведки был представлен в монографии [5]. Метод SVD успешно используется при анализе данных сейсморазведки [17, 6], гравиразведки [11], электроразведки [10] и комплексной интерпретации геофизических материалов [2].

Рисунок 1. Схема SVD квадратной матрицы коэффициентов СЛАУ при r < k

Выполнив замену переменных z = V^Tm, при умножении СЛАУ (1) слева на матрицу U^T получаем эквивалентную систему уравнений:

Sz = d, d = U^TV_R. (3)

Чтобы обеспечить устойчивое решение СЛАУ, достаточно вычислять только $r<k$ диагональных элементов матрицы S, соответствующих числу наибольших сингулярных чисел $δ_r>δ_{min}$ , где $δ_{min}$ — выбранное пороговое значение. Таким образом реализуется т. н. усеченный метод SVD (Truncated Singular Values Decomposition, TSVD), широко использующийся при сжатии изображений [13]. Обоснованный выбор параметра $γ$ в TSVD обеспечивает устойчивое приближенное регуляризованное решение СЛАУ (3), однако формализовать этот выбор весьма затруднительно.

Можно также использовать другой, редко применяющийся вариант регуляризации SVD, предложенный в работе [4]. Суть его заключается в переходе от СЛАУ (3) к нормальной системе уравнений, которая всегда разрешима:

S^TSz = S^Td. (4)

Для нее применима регуляризация в виде:

S^TSz + αz = S^Td. (5)

что позволяет выбирать параметр регуляризации α известными способами, в зависимости от априорной информации о погрешности исходных данных [19]. В случае работы со спутниковыми данными можно опираться на величину невязки ‖V_R - Gm‖_L² ≅ (1 - 2) мГал, при использовании материалов аэрогравиметрических или наземных точность регуляризованного решения СЛАУ (5) должна быть сопоставима со среднеквадратической погрешностью определения аномалий в редукции Буге.

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

При построении аналитической модели поля важную роль играет пространственное размещение масс: установлено, что разные варианты расположения эквивалентных источников, одинаково хорошо аппроксимирующие аномалии силы тяжести на поверхности Земли, могут генерировать заметно различающиеся между собой значения трансформант. Это объясняется использованием интеграла Пуассона для решения внешней задачи Дирихле при конечных пределах области интегрирования D и дискретном характере задания значений поля. В данном случае осуществляется совместная аппроксимация двух массивов данных, отвечающих непосредственно изучаемой площади D₁ и обрамляющей ее территории области D₂. Под трансформантой функции $u(φ,λ,R)$ подразумевается новая функция $v(φ,λ,R)$ , определенная на произвольной поверхности S в пределах области D₁: $v(φ,λ,R)=T\{u(φ,λ,R)\}=T\{u_1 (φ,λ,R)\}+T\{u_2 (φ,λ,R)\}$ , где Т — некоторый линейный оператор трансформации, а функция $u$ представляет собой сумму двух составляющих: $u=u_1+u_2$ , где $u_1$ — низкочастотная составляющая, моделируемая источниками нижнего уровня, расположенными в пределах области D₂, а $u_2$ — высокочастотная составляющая, созданная набором источников верхнего уровня, находящимися внутри области исследований D₁ [9].

Для расчетов использованы данные глобальной модели гравитационного поля Земли (ГПЗ) EIGEN-GRGS.RL04.MEAN-FIELD, полученной на основе данных спутниковых миссий GRACE и SLR в 2019 г. Источником информации о рельефе земной поверхности являлась модель ETOPO1. Диапазон изменения аномалий силы тяжести ~330 мГал, высотные отметки H лежат в пределах от −18 до 2418 м. Область D₁ размером ~5.15 млн. кв. км охватывает всю Сибирскую платформу и структуры ее обрамления (рис. 2). Шаг сети — 0.4°, размерность GRID-модели — 51 × 121 (N = 6171 точек). Эквивалентные источники располагаются на 45 км ниже уровня моря. Поле в области D₂ с координатами 48°–76° с. ш., 80°–136° в. д. задано с шагом 1°, глубина точечных масс составляет 111.3 км (N = 1653 точки).

Выбор относительных глубин источников, отвечающих шагу сети задания поля по меридиану для двух цифровых моделей поля V_R, влечет за собой сравнительно высокие числа обусловленности СЛАУ: cond G = 248012 при определении масс нижнего уровня и cond G = 9645 при определении масс верхнего уровня (обе матрицы являются невырожденными).

Рисунок 2. Карты изолиний аномалий силы тяжести (а) и высот рельефа дневной поверхности (б) Сибирской платформы

Результаты решения СЛАУ вида (1) в процессе двухэтапного процесса подбора масс m представлены в таблицах 1, 2. Относительные затраты времени на выполнение процедур матричной декомпозиции показаны на рисунке 3.

Таблица 1. Результаты решения СЛАУ для масс нижнего уровня (размер G — 1653 × 7824)

Тип разложения	Показатели качества решения СЛАУ, мГал				Время, с	Норма решения ‖m‖
Тип разложения	Минимум	Максимум	Среднее	СКО	Время, с	Норма решения ‖m‖
Холецкого	−84.73	86.66	0.018	±20.15	0.37	3970026
QR	−84.73	86.66	0.018	±20.15	1.75	3970026
SDV	−84.73	86.66	0.018	±20.18	10.91	1437685

Таблица 2. Результаты решения СЛАУ для масс верхнего уровня (размер G — 6171 × 6171)

Тип разложения	Показатели качества решения СЛАУ, мГал				Время, с	Норма решения ‖m‖
Тип разложения	Минимум	Максимум	Среднее	СКО	Время, с	Норма решения ‖m‖
Холецкого	−4.4×10^-7	4.9×10^-7	1.8×10^-14	3.4×10^-10	9.54	1979794
QR	−1.7×10^-11	2.4×10^-11	9.6×10^-15	1.1×10^-12	15.66	1979794
SDV	−1.7×10^-11	8.5×10^-12	4.8×10^-14	1.3×10^-12	74.74	1979595

Рисунок 3. Относительные затраты времени в процентах на выполнение процедур матричной декомпозиции для нижних (а) и верхних (б) масс: 1 — разложение Холецкого; 2 — QR-разложение; 3 — SVD

Как очевидно, все виды матричной декомпозиции обеспечивают высокое качество построения аналитической модели поля при заметных различиях в скорости вычислений. Нужно пояснить, что при определении масс нижнего уровня методом SVD использовалось 95 % сингулярных чисел матрицы S и уменьшение ‖m‖ обусловлено применением этой регуляризации (рис. 4).

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

Пересчет гравитационного поля на фрагмент внешней по отношению к Земле сферической поверхности с радиусом $R'= R_0+H$ , где H — нормальная высота, в данном случае отвечает приведению поля $V_Z$ на горизонтальную плоскость $z=const$ для плоской Земли. При работе в сферических координатах $∑φ,λ,r$ используются достаточно простые аналитические выражения для вторых производных гравитационного потенциала V [7]. Аналогом вертикальной производной силы тяжести $V_{ZZ}$ является 2-я радиальная производная силы тяжести $V_{RR}$ . Соответствующую аналогию можно также провести между вторыми производными $V_{ZX}$ , $V_{ZY}$ и $V_{Rλ}$ , $V_{Rφ}$ . На рисунке 5 приведены карта региональной составляющей гравитационного поля $V_R$ изучаемой территории, полученная путем пересчета в верхнее полупространство на высоту 250 км, и карта локальной (остаточной) составляющей этого поля.

Рисунок 5. Карты изоаномал региональной (а) и локальной (б) составляющих гравитационного поля Сибирской платформы

Объемное распределение эффективного параметра, пропорционального плотности горных пород, построено в системе КОСКАД-3D с применением статистических, спектрально-корреляционных методов и алгоритма адаптивной фильтрации в окне живой формы на основе преобразования Б. А. Андреева [16]. Блок-диаграмма квазиплотности (рис. 6) представляет собой результат перехода от 2D- к 3D-геоизображениям на основе аналитических компьютерных процедур логико-математической генерализации. Последовательное повышение уровня генерализации обеспечивает проявление на геоизображениях элементов все более крупных геосистем, что способствует выявлению качественно новой информации и выявлению скрытых закономерностей [3]. Построенная блок-диаграмма характеризует основные особенности глубинного строения Сибирской платформы и обрамляющих ее структур, в т. ч. слабо проявленные на карте изоаномал наблюденного гравитационного поля.

Рисунок 6. Блок-диаграмма относительного распределения плотности в литосфере Сибирской платформы

Выводы

Рассмотренные способы матричной декомпозиции могут успешно использоваться при построении аналитических моделей региональных аномалий силы тяжести и их трансформации. В первую очередь это относится к Арктике и Антарктике, где уменьшение длины параллелей влечет за собой плохую обусловленность СЛАУ, решение которых требуется для определения эквивалентных масс [8]. SVD может рассматриваться как альтернатива классическим методам регуляризации, при числе точек поля $m=n≥200$ этот подход обладает заметными преимуществами по скорости вычислений [15], хотя и является самым медленным из рассмотренных. Наибольшую скорость вычислений обеспечивает разложение Холецкого, но при очень малых (близких к машинному нулю) коэффициентах G системы уравнений эта процедура может отказать или потерять устойчивость. QR-разложение является менее быстрым, но более надежным способом решения плохо обусловленных нормальных систем уравнений, т. к. сохраняет евклидову норму матрицы коэффициентов. Рассмотренные в докладе алгоритмы ориентированы на извлечение геологической информации, представленной в спутниковых (альтиметрических) глобальных моделях ГПЗ, но без ограничений могут работать с другими источниками данных.

Благодарности

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

1. Бахвалов Н. С. Численные методы / Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. — М.: Наука, 2000.

2. Белашев Б. З. Моделирование литосферы региона Белого моря с использованием декомпозиции аномальных гравитационного и магнитного полей / Белашев Б. З., Бакунович Л. И., Шаров Н. В. // Геодинамика и тектонофизика. — 2023. — Т. 14. — № 5:0720.

3. Берлянт А. М. Теория геоизображений / Берлянт А. М. — М.: ГЕОС, 2006.

4. Васин В. В. Основы теории некорректных задач / Васин В. В. — Новосибирск: изд-во СО РАН, 2020.

5. Вахромеев Г. С. Моделирование в разведочной геофизике / Вахромеев Г. С., Давыденко А. Ю. — М.: Недра, 1987.

6. Гадыльшин К. Г. Обращение полных волновых полей нелинейным методом наименьших квадратов: SVD анализ / Гадыльшин К. Г., Чеверда В. А. // Вычислительные методы и программирование. — 2014. — Т. 15. — № 3. — С. 499–513.

7. Долгаль А. С. Моделирование аномалий силы тяжести системой точечных масс на сфероообразной Земле / Долгаль А. С., Костицын В. И., Новикова [и др.] // Геофизика. — 2023. — № 5. — С. 10–17.

8. Долгаль А. С. Моделирование аномалий силы тяжести эквивалентными источниками в полярных областях Земли / Долгаль А. С., Огородова И. В., Рыжов Н. В., Христенко Л. А. // Вестник Пермского университета. Геология. — 2025. — Т. 24. — № 3. — С. 226–235.

9. Долгаль А. С. Повышение точности трансформации региональных гравитационных аномалий: учет влияния масс, находящихся за пределами изучаемой площади / Долгаль А. С., Рыжов Н. В., Хохлова В. В. // Вестник АГГИ. — 2025. — № 2(2). — С. 23–34.

10. Кабанихин С. И. Сингулярное разложение в задаче об источнике / Кабанихин С. И., Криворотько О. И. // Сибирский журнал вычислительной математики. — 2012. — Т. 15. — № 2. — С. 205–211.

11. Канайкин В. С. Особенности трансформации и инверсии гравитационного поля с применением дисперсионного и регрессионного анализов / Канайкин В. С., Турутанов Е. X., Буянтогтох Б. // Известия Сибирского отделения секции наук о Земле Российской академии естественных наук. Геология, разведка и разработка месторождений полезных ископаемых. — 2018. — Т. 41. — № 3(64). — С. 93–105.

12. Карпик А. П. Исследования спектральных характеристик глобальных моделей гравитационного поля Земли, полученных по космическим миссиям CHAMP, GRAСE и GOCE / Карпик А. П., Канушин В. Ф., Ганагина И. Г., Голдобин Д. Н., Мазурова Е. М. // Гироскопия и навигация. — 2014. — № 4(87). — С. 34–44.

13. Каханер Д. Численные методы и программное обеспечение (пер. с англ. под ред. Х. Д. Икрамова). — М.: Мир, 2001.

14. Конешов В. Н. Особенности сравнительной оценки глобальных моделей гравитационного поля Земли / Конешов В. Н., Непоклонов В. Б., Спиридонова Е. С., Максимова М. В. // Физика Земли. — 2020. — № 2. — С. 115–126.

15. Леонов А. С. Решение некорректно поставленных обратных задач: очерк теории, практические алгоритмы и демонстрации в МАТЛАБ / Леонов А. С. — М.: ЛЕНАНД, 2024.

16. Петров А. В. Компьютерная технология статистического и спектрально-корреляционного анализа данных КОСКАД 3D и практические результаты / Петров А. В., Демура Г. В., Зиновкин С. В. // Недропользование XXI век. — 2017. — № 1(64). — С. 44–59.

17. Сильвестров И.Ю. Анализ сингулярного разложения линеаризованного оператора динамической теории упругости для случая вертикального сейсмического профилирования / Сильвестров И. Ю. // Вычислительные технологии. — 2007. — Т. 12. — № 6. — С. 90–100.

18. Страхов В. Н. Разработка теории и компьютерной технологии построения линейных аналитических аппроксимаций гравитационных и магнитных полей / Страхов В. Н., Керимов И. А., Степанова И. Э. — М.: ОИФЗ РАН, 2009.

19. Тихонов А. Н. Методы решения некорректных задач: учебное пособие. — М.: ЛЕНАНД, 2022.

Отправить рукопись Скачать PDF
Текст JATS XML