ВведениеПод аналитической моделью аномалий силы тяжести ∆g∆g в данном случае подразумевается множество параметров сеточной модели эквивалентных источников (точечных масс), порождающих модельное гравитационное поле ∆gT∆g_T, весьма близкое наблюденному полю ‖∆g-∆gT‖≤ε‖∆g-∆g_T ‖≤ε, где величина εε определяется предполагаемой нормой помех δδ, δ=‖ε‖δ=‖ε‖. В это множество входят геометрические и физические параметры, т. е. априорно заданные координаты размещения источников и значения их масс соответственно. Будем рассматривать моделирование аномалий силы тяжести в редукции Буге на фрагменте сферообразной Земли. По классификации В. Н. Страхова, речь пойдет о региональном варианте линейных аналитических аппроксимаций, который применяется для сравнительно больших территорий размером 105–106 км2 и более [18]. В этом случае в качестве модели Земли используется сфера Каврайского, для расчетов выполняется переход от геодезических координат (широты B и долготы L) к сферическим координатам ∑φλR∑R , а аномалии силы тяжести ∆g∆g отождествляются с 1-й радиальной производной гравитационного потенциала VRV_R [7].Значения аппроксимирующих масс m определяются путем решения системы линейных алгебраических уравнений (СЛАУ):Gm = VR,       (1)где G — квадратная матрица значений гравитационных эффектов для точечного источника с единичной массой (m=1)(m=1), VR — вектор наблюденных аномалий силы тяжести.Если при решении практических задач порядок СЛАУ (1) составляет N ≌ (1 – 5) ´ 104 и менее, то для ее решения можно отказаться от использования приближенных итерационных методов [1] и перейти к алгоритмам, основанным на матричной декомпозиции [13]. Этот вариант достаточно часто возникает при работе с глобальными моделями гравитационного поля Земли [12, 14], разрешающая способность которых ограничивает число независимых значений аномалий силы тяжести в пределах изучаемых территорий. Далее будет рассмотрено применение разложения Холецкого, QR-разложения и сингулярного разложения (SVD — singular values decomposition) в процессе истокообразной аппроксимации региональных гравитационных аномалий VRV_R в редукции Буге. Суть данного подхода заключается в представлении матрицы G в виде произведения двух и более матриц, обладающих  определенными свойствами (ортогональностью, симметричностью, диагональностью), облегчающими решение задачи (1).Методы разложения матрицМетоды матричной декомпозиции, по мнению, высказанному в журнале SIAM Rewiev в 2000 г., входят в список из 10 алгоритмов, оказавших наибольшее влияние на науку и индустрию в ХХ веке. Фундаментальные работы по этой проблеме были опубликованы учеными из США Полом Дуайером (Paul Dwyer) и Альстоном Скоттом Хаусхолдером (Alston Scott Householder) в 1951 и 1954 гг. Далее рассмотрим 3 способа факторизации матриц: разложение Холецкого, QR-разложение и сингулярное разложение.Разложение Холецкого позволяет представить положительно определенную симметричную матрицу коэффициентов нормальной системы уравнений в виде произведения левой и правой треугольных матриц GTG = LLT, где L — нижняя треугольная матрица со строго положительными элементами главной диагонали. QR-разложение — это представление невырожденной матрицы G в виде G = QR, где Q — ортогональная матрица (QT = Q-1), a R — верхняя треугольная матрица [13]. QR-разложение выполняется с сохранением нормы Фробениуса исходной матрицы. Это является его преимуществом по сравнению с классическим подходом в методе наименьших квадратов, где требуется решение нормальной СЛАУ вида GTGm = GTVR. Число обусловленности матрицы коэффициентов такой системы достаточно велико: cond GTG ≈ cond2 G .Сингулярное разложение (Singular Values Decomposition, SVD) позволяет представить матрицу G размером k×kk×k в виде, представленном на рисунке 1.G = USVT,       (2)где U, V — ортогональные матрицы (UT = U-1, VT = V-1), (U^T=U^{-1},V^T=V^{-1})S ={δii}=\{δ_{ii} \} — диагональная матрица с элементами δi(δi≥δi+1)δ_i (δ_i≥δ_{i+1}), которые называются сингулярными числами матрицы [13].Первый в нашей стране опыт применения сингулярного разложения для факторного анализа петрофизических данных и решения линейной обратной задачи гравиразведки был представлен в монографии [5]. Метод SVD успешно используется при анализе данных сейсморазведки [17, 6], гравиразведки [11], электроразведки [10] и комплексной интерпретации геофизических материалов [2].Рисунок 1. Схема SVD квадратной матрицы коэффициентов СЛАУ при r < kВыполнив замену переменных z = VTm, при умножении СЛАУ (1) слева на матрицу UT получаем эквивалентную систему уравнений:Sz = d, d = UTVR.       (3)Чтобы обеспечить устойчивое решение СЛАУ, достаточно вычислять только r<kr<k диагональных элементов матрицы S, соответствующих числу наибольших сингулярных чисел δr>δminδ_r>δ_{min}, где δminδ_{min} — выбранное пороговое значение. Таким образом реализуется т. н. усеченный метод SVD (Truncated Singular Values Decomposition, TSVD), широко использующийся при сжатии изображений [13]. Обоснованный выбор параметра γγ в TSVD обеспечивает устойчивое приближенное регуляризованное решение СЛАУ (3), однако формализовать этот выбор весьма затруднительно.Можно также использовать другой, редко применяющийся вариант регуляризации SVD, предложенный в работе [4]. Суть его заключается в переходе от СЛАУ (3) к нормальной системе уравнений, которая всегда разрешима:STSz = STd.       (4)Для нее применима регуляризация в виде:STSz + αz = STd.       (5)что позволяет выбирать параметр регуляризации α известными способами, в зависимости от априорной информации о погрешности исходных данных [19]. В случае работы со спутниковыми данными можно опираться на величину невязки ‖VR - Gm‖L2 ≅ (1 - 2) мГал, при использовании материалов аэрогравиметрических или наземных точность регуляризованного решения СЛАУ (5) должна быть сопоставима со среднеквадратической погрешностью определения аномалий в редукции Буге.Построение аналитической модели гравитационного поля Сибирской платформыПри построении аналитической модели поля важную роль играет пространственное размещение масс: установлено, что разные варианты расположения эквивалентных источников, одинаково хорошо аппроксимирующие аномалии силы тяжести на поверхности Земли, могут генерировать заметно различающиеся между собой значения трансформант. Это объясняется использованием интеграла Пуассона для решения внешней задачи Дирихле при конечных пределах области интегрирования D и дискретном характере задания значений поля. В данном случае осуществляется совместная аппроксимация двух массивов данных, отвечающих непосредственно изучаемой площади D1 и обрамляющей ее территории области D2. Под трансформантой функции u(φ,λ,R)u(φ,λ,R) подразумевается новая функция v(φ,λ,R)v(φ,λ,R), определенная на произвольной поверхности S в пределах области D1: v(φ,λ,R)=T{u(φ,λ,R)}=T{u1(φ,λ,R)}+T{u2(φ,λ,R)}v(φ,λ,R)=T\{u(φ,λ,R)\}=T\{u_1 (φ,λ,R)\}+T\{u_2 (φ,λ,R)\}, где Т — некоторый линейный оператор трансформации, а функция uu представляет собой сумму двух составляющих: u=u1+u2u=u_1+u_2, где u1u_1 — низкочастотная составляющая, моделируемая источниками нижнего уровня, расположенными в пределах области D2, а u2u_2 — высокочастотная составляющая, созданная набором источников верхнего уровня, находящимися внутри области исследований D1 [9].Для расчетов использованы данные глобальной модели гравитационного поля Земли (ГПЗ) EIGEN-GRGS.RL04.MEAN-FIELD, полученной на основе данных спутниковых миссий GRACE и SLR в 2019 г. Источником информации о рельефе земной поверхности являлась модель ETOPO1. Диапазон изменения аномалий силы тяжести ~330 мГал, высотные отметки H лежат в пределах от −18 до 2418 м. Область D1 размером ~5.15 млн. кв. км охватывает всю Сибирскую платформу и структуры ее обрамления (рис. 2). Шаг сети — 0.4°, размерность GRID-модели — 51 × 121 (N = 6171 точек). Эквивалентные источники располагаются на 45 км ниже уровня моря. Поле в области D2 с координатами 48°–76° с. ш., 80°–136° в. д. задано с шагом 1°, глубина точечных масс составляет 111.3 км (N = 1653 точки).Выбор относительных глубин источников, отвечающих шагу сети задания поля по меридиану для двух цифровых моделей поля VR, влечет за собой сравнительно высокие числа обусловленности СЛАУ: cond G = 248012 при определении масс нижнего уровня и cond G = 9645 при определении масс верхнего уровня (обе матрицы являются невырожденными).Рисунок 2. Карты изолиний аномалий силы тяжести (а) и высот рельефа дневной поверхности (б) Сибирской платформыРезультаты решения СЛАУ вида (1) в процессе двухэтапного процесса подбора масс m представлены в таблицах 1, 2. Относительные затраты времени на выполнение процедур матричной декомпозиции показаны на рисунке 3.Таблица 1. Результаты решения СЛАУ для масс нижнего уровня (размер G — 1653 × 7824)Тип разложенияПоказатели качества решения СЛАУ, мГалВремя, сНорма решения‖m‖МинимумМаксимумСреднееСКОХолецкого−84.7386.660.018±20.150.373970026QR−84.7386.660.018±20.151.753970026SDV−84.7386.660.018±20.1810.911437685 Таблица 2. Результаты решения СЛАУ для масс верхнего уровня (размер G — 6171 × 6171)Тип разложенияПоказатели качества решения СЛАУ, мГалВремя, сНорма решения ‖m‖МинимумМаксимумСреднееСКОХолецкого−4.4×10-74.9×10-71.8×10-143.4×10-109.541979794QR−1.7×10-112.4×10-119.6×10-151.1×10-1215.661979794SDV−1.7×10-118.5×10-124.8×10-141.3×10-1274.741979595 Рисунок 3. Относительные затраты времени в процентах на выполнение процедур матричной декомпозиции для нижних (а) и верхних (б) масс: 1 — разложение Холецкого; 2 — QR-разложение; 3 — SVDКак очевидно, все виды матричной декомпозиции обеспечивают высокое качество построения аналитической модели поля при заметных различиях в скорости вычислений. Нужно пояснить, что при определении масс нижнего уровня методом SVD использовалось 95 % сингулярных чисел матрицы S и уменьшение ‖m‖ обусловлено применением этой регуляризации (рис. 4).Рисунок 4. Сингулярные числа δ матрицы коэффициентов: моделирование аномалий силы тяжести Сибирской платформы с усечением S при r = 1488Пересчет гравитационного поля на фрагмент внешней по отношению к Земле сферической поверхности с радиусом R′=R0+HR&#39;= R_0+H, где H — нормальная высота, в данном случае отвечает приведению поля VZV_Z на горизонтальную плоскость z=constz=const для плоской Земли. При работе в сферических координатах ∑φ,λ,r∑φ,λ,r используются достаточно простые аналитические выражения для вторых производных гравитационного потенциала V [7]. Аналогом вертикальной производной силы тяжести VZZV_{ZZ} является 2-я радиальная производная силы тяжести VRRV_{RR}. Соответствующую аналогию можно также провести между вторыми производными VZXV_{ZX}, VZYV_{ZY} и VRλV_{Rλ}, VRφV_{Rφ}. На рисунке 5 приведены карта региональной составляющей гравитационного поля VRV_R изучаемой территории, полученная путем пересчета в верхнее полупространство на высоту 250 км, и карта локальной (остаточной) составляющей этого поля.Рисунок 5. Карты изоаномал региональной (а) и локальной (б) составляющих гравитационного поля Сибирской платформыОбъемное распределение эффективного параметра, пропорционального плотности горных пород, построено в системе КОСКАД-3D с применением статистических, спектрально-корреляционных методов и алгоритма адаптивной фильтрации в окне живой формы на основе преобразования Б. А. Андреева [16]. Блок-диаграмма квазиплотности (рис. 6) представляет собой результат перехода от 2D- к 3D-геоизображениям на основе аналитических компьютерных процедур логико-математической генерализации. Последовательное повышение уровня генерализации обеспечивает проявление на геоизображениях элементов все более крупных геосистем, что способствует выявлению качественно новой информации и выявлению скрытых закономерностей [3]. Построенная блок-диаграмма характеризует основные особенности глубинного строения Сибирской платформы и обрамляющих ее структур, в т. ч. слабо проявленные на карте изоаномал наблюденного гравитационного поля.Рисунок 6. Блок-диаграмма относительного распределения плотности в литосфере Сибирской платформыВыводыРассмотренные способы матричной декомпозиции могут успешно использоваться при построении аналитических моделей региональных аномалий силы тяжести и их трансформации. В первую очередь это относится к Арктике и Антарктике, где уменьшение длины параллелей влечет за собой плохую обусловленность СЛАУ, решение которых требуется для определения эквивалентных масс [8]. SVD может рассматриваться как альтернатива классическим методам регуляризации, при числе точек поля m=n≥200m=n≥200  этот подход обладает заметными преимуществами по скорости вычислений [15], хотя и является самым медленным из рассмотренных. Наибольшую скорость вычислений обеспечивает разложение Холецкого, но при очень малых (близких к машинному нулю) коэффициентах G системы уравнений эта процедура может отказать или потерять устойчивость. QR-разложение является менее быстрым, но более надежным способом решения плохо обусловленных нормальных систем уравнений, т. к. сохраняет евклидову норму матрицы коэффициентов. Рассмотренные в докладе алгоритмы ориентированы на извлечение геологической информации, представленной в спутниковых (альтиметрических) глобальных моделях ГПЗ, но без ограничений могут работать с другими источниками данных.БлагодарностиИсследование выполнено при финансовой поддержке Министерства науки и высшего образования РФ в рамках государственного задания (рег. номер НИОКТР: 126012716041-5).