Application of matrix decomposition to construct analytical models  of gravity anomalies

Aleksandr Dolgal&#039;

doi:doi:

Home / Journals / Vestnik AGGE / Application of matrix decomposition to construct analytical models of gravity anomalies

Submit manuscript Download PDF
Text

To cite

Citations:

Application of matrix decomposition to construct analytical models of gravity anomalies

Journal: VESTNIK AGGE

Rubrics: ENGINEERING AND MINING GEOPHYSICS

Aleksandr Dolgal' ¹

Author and publication information

Authors:

1. Mining Institute of the Ural Branch of the Russian Academy of Sciences (EGPGlavnyy nauchnyy sotrudnik)
employee from 01.01.2002 until now

Ural State Mining University

Ekaterinburg, Russian Federation

Type:

Article

Pages:

from 1 to 1

Status:

Subject area:

VAK Russia 1.6
VAK Russia 1.2.1
VAK Russia 1.2.2
UDC 55
UDC 550
UDC 550.3

Language:

Russian

Keywords:

gravity field, Bouguer anomalies, analytical field model, point masses, system of linear equations, Cholesky decomposition, singular value decomposition, QR decomposition, regularization, Arctic, Antarctica.

Funding:

The study was carried out with the financial support of the Ministry of Science and Higher Education of the Russian Federation within the framework of a state assignment (reg. number of R&D: 126012716041-5).

Abstract and keywords

Abstract:
The application of matrix decomposition to constructing analytical models of regional gravity anomalies using point masses on a spherical Earth is considered. The point masses are placed at two depth levels. The Cholesky, SVD, and QR decompositions are effectively used to solve systems of linear algebraic equations in the process of sourcewise approximation. An example of modeling the gravity field in the Bouguer reduction is presented for the territory of the Siberian Platform, approximately 5 million square kilometers in size, with coordinates 52°–72° N, 84°–132° E.

Keywords:
gravity field, Bouguer anomalies, analytical field model, point masses, system of linear equations, Cholesky decomposition, singular value decomposition, QR decomposition, regularization, Arctic, Antarctica.

Text

Text (PDF): Read Download

Введение

Под аналитической моделью аномалий силы тяжести $∆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-разложение является менее быстрым, но более надежным способом решения плохо обусловленных нормальных систем уравнений, т. к. сохраняет евклидову норму матрицы коэффициентов. Рассмотренные в докладе алгоритмы ориентированы на извлечение геологической информации, представленной в спутниковых (альтиметрических) глобальных моделях ГПЗ, но без ограничений могут работать с другими источниками данных.

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

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

References

1. Bahvalov N. S. Chislennye metody / Bahvalov N. S., Zhidkov N. P., Kobel'kov G. M. — M.: Nauka, 2000 (in Russ.).

2. Belashev B. Z. Modelirovanie litosfery regiona Belogo morya s ispol'zovaniem dekompozicii anomal'nyh gravitacionnogo i magnitnogo poley / Belashev B. Z., Bakunovich L. I., Sharov N. V. // Geodinamika i tektonofizika. — 2023. — Vol. 14. — № 5:0720 (in Russ.).

3. Berlyant A. M. Teoriya geoizobrazheniy / Berlyant A. M. — M.: GEOS, 2006 (in Russ.).

4. Vasin V. V. Osnovy teorii nekorrektnyh zadach / Vasin V. V. — Novosibirsk: izd-vo SO RAN, 2020 (in Russ.).

5. Vahromeev G. S. Modelirovanie v razvedochnoy geofizike / Vahromeev G. S., Davydenko A. Yu. — M.: Nedra, 1987 (in Russ.).

6. Gadyl'shin K. G. Obraschenie polnyh volnovyh poley nelineynym metodom naimen'shih kvadratov: SVD analiz / Gadyl'shin K. G., Cheverda V. A. // Vychislitel'nye metody i programmirovanie. — 2014. — Vol. 15. — № 3. — P. 499–513 (in Russ.).

7. Dolgal' A. S. Modelirovanie anomaliy sily tyazhesti sistemoy tochechnyh mass na sferooobraznoy Zemle / Dolgal' A. S., Kosticyn V. I., Novikova [et al.] // Geofizika. — 2023. — № 5. — P. 10–17 (in Russ.).

8. Dolgal' A. S. Modelirovanie anomaliy sily tyazhesti ekvivalentnymi istochnikami v polyarnyh oblastyah Zemli / Dolgal' A. S., Ogorodova I. V., Ryzhov N. V., Hristenko L. A. // Vestnik Permskogo universiteta. Geologiya. — 2025. — Vol. 24. — № 3. — P. 226–235 (in Russ.).

9. Dolgal' A. S. Povyshenie tochnosti transformacii regional'nyh gravitacionnyh anomaliy: uchet vliyaniya mass, nahodyaschihsya za predelami izuchaemoy ploschadi / Dolgal' A. S., Ryzhov N. V., Hohlova V. V. // Vestnik AGGI. — 2025. — № 2(2). — P. 23–34 (in Russ.).

10. Kabanihin S. I. Singulyarnoe razlozhenie v zadache ob istochnike / Kabanihin S. I., Krivorot'ko O. I. // Sibirskiy zhurnal vychislitel'noy matematiki. — 2012. — Vol. 15. — № 2. — P. 205–211 (in Russ.).

11. Kanaykin V. S. Osobennosti transformacii i inversii gravitacionnogo polya s primeneniem dispersionnogo i regressionnogo analizov / Kanaykin V. S., Turutanov E. X., Buyantogtoh B. // Izvestiya Sibirskogo otdeleniya sekcii nauk o Zemle Rossiyskoy akademii estestvennyh nauk. Geologiya, razvedka i razrabotka mestorozhdeniy poleznyh iskopaemyh. — 2018. — Vol. 41. — № 3(64). — P. 93–105 (in Russ.).

12. Karpik A. P. Issledovaniya spektral'nyh harakteristik global'nyh modeley gravitacionnogo polya Zemli, poluchennyh po kosmicheskim missiyam CHAMP, GRASE i GOCE / Karpik A. P., Kanushin V. F., Ganagina I. G., Goldobin D. N., Mazurova E. M. // Giroskopiya i navigaciya. — 2014. — № 4(87). — P. 34–44 (in Russ.).

13. Kahaner D. Chislennye metody i programmnoe obespechenie. — M.: Mir, 2001 (in Russ.).

14. Koneshov V. N. Osobennosti sravnitel'noy ocenki global'nyh modeley gravitacionnogo polya Zemli / Koneshov V. N., Nepoklonov V. B., Spiridonova E. S., Maksimova M. V. // Fizika Zemli. — 2020. — № 2. — P. 115–126 (in Russ.).

15. Leonov A. S. Reshenie nekorrektno postavlennyh obratnyh zadach: ocherk teorii, prakticheskie algoritmy i demonstracii v MATLAB / Leonov A. S. — M.: LENAND, 2024 (in Russ.).

16. Petrov A. V. Komp'yuternaya tehnologiya statisticheskogo i spektral'no-korrelyacionnogo analiza dannyh KOSKAD 3D i prakticheskie rezul'taty / Petrov A. V., Demura G. V., Zinovkin S. V. // Nedropol'zovanie XXI vek. — 2017. — № 1(64). — P. 44–59 (in Russ.).

17. Sil'vestrov I.Yu. Analiz singulyarnogo razlozheniya linearizovannogo operatora dinamicheskoy teorii uprugosti dlya sluchaya vertikal'nogo seysmicheskogo profilirovaniya / Sil'vestrov I. Yu. // Vychislitel'nye tehnologii. — 2007. — Vol. 12. — № 6. — P. 90–100 (in Russ.).

18. Strahov V. N. Razrabotka teorii i komp'yuternoy tehnologii postroeniya lineynyh analiticheskih approksimaciy gravitacionnyh i magnitnyh poley / Strahov V. N., Kerimov I. A., Stepanova I. E. — M.: OIFZ RAN, 2009 (in Russ.).

19. Tihonov A. N. Metody resheniya nekorrektnyh zadach: uchebnoe posobie. — M.: LENAND, 2022 (in Russ.).

Submit manuscript Download PDF
Text JATS XML