ВведениеВ прикладной геофизике при обработке и интерпретации материалов гравиметрической съемки обычно используется модель «плоской Земли» и отвечающая ей декартова система прямоугольных координат. Значительно реже находит применение при изучении больших территорий модель «сферической Земли» [1]. Некоторые исследователи работают с более сложными моделями «эллипсоидальной Земли», являющимися более точным приближением формы нашей планеты [2].Рассмотрим экспериментальные оценки различий результатов пересчета аномалий силы тяжести в верхнее полупространство для перечисленных выше моделей. Пересчет на высоту 50 км выполнялся путем подбора масс эквивалентных источников с последующим решением прямой задачи гравиразведки от полученной аппроксимационной конструкции [3]. Эквивалентными источниками во всех случаях являлись шары (точечные массы), находящиеся под точками задания поля на фиксированной глубине h0, близкой к шагу сети цифровой GRID-модели.При сопоставлении «плоской» и «сферической» моделей Земли исходными данными являлись цифровые GRID-модели гравитационного поля ΔΔg в редукции Буге и нормальных высот рельефа земной поверхности в пределах листа О-40 карты 1:1 000 000 масштаба (Пермский край). Это предельный размер территории (т. н. сферического двуугольника), при котором целесообразно использовать систему плоских прямоугольных координат Гаусса — Крюгера и соответствующую ей модель Земли. Сеть точек составила 8 ×× 8 км, размер матриц — 63 строки, 53 столбца. Глубина источников h0 = 9.6 км, при этом была достигнута точность решения СЛАУ F2 = 0.001 мГал.Плановое расположение точек поля было одинаковым для обеих моделей Земли, но в первом случае уровенной поверхностью являлась горизонтальная плоскость z = 0, во втором — фрагмент поверхности сферы, обладающей радиусом R = 6371.1 км. Две этих уровенных поверхности соприкасались между собой в юго-западном углу территории. Шар радиуса 6371.1 км по своим размерам, площади поверхности и объему близок к земному эллипсоиду. Для «плоской» Земли аномалии силы тяжести отвечают вертикальной производной VZ гравитационного потенциала V, которую обычно принято обозначать ΔΔg. Для модели «сферической» Земли используется уже радиальная производная потенциала VR. Соответственно, пересчет поля выполнялся на плоскость z = 50 км и на сферу радиусом R = 6421.1 км, а также определялись значения разностного поля δ = VR - VZ (таб. 1, рис. 1).Разностное поле характеризуется среднеквадратическим отклонением (СКО) ~~ 3.4 % от СКО полей на высоте 50 км, диапазон его изменения составляет ~~ 4.6 % от размаха каждого из этих полей. Таким образом, можно говорить о значимых различиях результатов трансформации аномалий силы тяжести, полученных в рамках «плоской» и «сферической» моделей Земли. Степень отличия модели «сферической Земли» от эллипсоида Красовского существенно ниже: СКО разностного поля δδ составляет ~~ 0.04 % от СКО аномалий силы тяжести на высоте 50 км, а диапазон изменения δδ не превышает ~~ 0.08 % [4].Более точное приближение к форме Земли обеспечивает сфера В. В. Каврайского с радиусом R = 6372.9 км. Напомним, что при использовании «сферических» моделей Земли долгота L остается неизменной, но происходит замена геодезической широты B на геоцентрическую широту Ф. Направление силы тяжести теперь отвечает радиус-вектору, а не нормали ll к поверхности эллипсоида. Для перехода от геодезической широты к геоцентрической используется формула: Ф=B-8′39″sin2BФ = B - 8' 39'' sin 2B [5]. Использование «квазиэллипсоидальной Земли» — сферы Каврайского (при замене Ф на кошироту ϕϕ) позволяет существенно ускорить вычислительный процесс, т. к. в этом случае не требуется перехода к топоцентрической системе координат в каждой точке задания поля [6].Таблица 1. Характеристики полей «плоской» и «сферической» моделей Земли на высоте 50 кмМодель ЗемлиСтатистические параметры, мГалМинимумМаксимумСреднееСКО«Плоская Земля»- 11.42614.3441.4225.875«Сферическая Земля»- 11.41814.2901.5945.854Разностное поле- 0.6260.5890.1200.199 Рисунок 1. Оценка различий моделей «плоской» и «сферической» Земли: а — гравитационное поле ΔΔg на земной поверхности; б — гравитационное поле VR на высоте 50 км; в — разностное поле δδ. Лист карты О-40 (Пермский край)Однако научно обоснованный выбор модели Земли не снимает всех сложностей, возникающих при трансформации региональных аномалий силы тяжести VR. Речь идет об учете влияния гравитационного эффекта масс, которые находятся за пределами территории исследований, а также о решении систем линейных алгебраических уравнений (СЛАУ) большой размерности с приближенно заданной правой частью. Далее в статье будут рассмотрены вопросы, связанные с этими сложностями.Учет влияния сторонних массИмеются многочисленные примеры, свидетельствующие о том, что разные пространственные распределения эквивалентных источников, одинаково хорошо аппроксимирующие аномалии силы тяжести на поверхности Земли, могут генерировать заметно различающиеся между собой значения трансформант [7]. С математической точки зрения это объясняется использованием интеграла Пуассона для решения внешней задачи Дирихле при конечных пределах области интегрирования и дискретном характере задания значений поля. В большинстве случаев наиболее заметные искажения трансформант возникают на периферии области задания поля D1, их можно назвать краевыми эффектами.В идеале избежать краевых эффектов можно при построении аппроксимационной конструкции, полностью адекватной реальному распределению аномалиеобразующих масс (плотности) в геологической среде, причем в объеме V2, существенно превышающем размер изучаемого фрагмента V1, V1 ⊂ V2 (рис. 2). На практике этого достичь невозможно, т. к., начиная с какого-то объема петрофизической информации, «лишаются смысла гравиразведочные (магниторазведочные) работы, цель которых состоит в изучении распределения плотности (намагниченности) в геологической среде» [8].Рисунок 2. Схема, поясняющая соотношение между объемом V1, содержащим аномалиеобразующие массы, и изучаемым фрагментом геоплотностной среды V2С целью минимизации краевых эффектов предлагается использовать два уровня глубин размещения эквивалентных источников (точечных масс): «приповерхностный» и «глубинный», а также информацию, относящуюся к самой области D1 и к существенно превышающей ее по площади сопредельной территории D2. Имеются примеры современных исследований, доказывающие возможность повышения точности расчета трансформант при размещении точечных масс на разных глубинах h от земной поверхности, в частности — работа [9]. Совокупность значений и пространственных координат этих масс представляет собой аналитическую модель интерпретируемого гравитационного поля VR (рис. 3). Высокая скорость вычислений обеспечивается за счет разрежения сети задания значений поля в пределах области D2 ⊃ D1.Рисунок 3. Схема расположения точечных масс в аналитической модели гравитационного поля: красный цвет — «приповерхностные» источники с глубиной h1, синий цвет — «глубинные» источники с глубиной h2 > h1Этот способ не является новым [10], но он использовался лишь для построения «локальных аппроксимаций» (термин В. Н. Страхова) в прямоугольной системе координат. В данном случае под линейной трансформантой функции u(φ,λ,R)u(φ, λ, R) подразумевается новая функция v(φ,λ,R)v(φ, λ, R), определенная на произвольной поверхности S в пределах области D1:v(φ,λ,R)=T{u(φ,λ,R)}=T{u1(φ,λ,R)}+T{u2(φ,λ,R)},v\left( \varphi,\lambda,R \right) = T\left\{ u\left( \varphi,\lambda,R \right) \right\} = T\left\{ u_{1}\left( \varphi,\lambda,R \right) \right\} + T\left\{ u_{2}\left( \varphi,\lambda,R \right) \right\},     (1)где Т — некоторый линейный оператор трансформации, а функция uu представляет собой сумму двух составляющих: u=u1+u2u = u_1 + u_2, где u1u_1 — модельное поле VR «глубинных» масс, а 𝑢2𝑢_2 — разность наблюденного и модельного полей внутри области исследований D1.Определение значений масс эквивалентных источниковМассы эквивалентных источников определяются с помощью решения системы линейных алгебраических уравнений (СЛАУ) с приближенно заданной правой частью вида:Gm=u,     (2)где G — квадратная матрица значений гравитационных эффектов элементарных источников с единичной массой (m = 1), m — вектор неизвестных значений аномальных масс, u — вектор значений аномального гравитационного поля VR в редукции Буге. Под гравитационным эффектом в данном случае подразумевается 1-я радиальная производная гравитационного потенциала точечного источника VR(R0,φ0,λ0)=f(R0-rcosω)r03V_{R}\left( R_{0},\varphi_{0},\lambda_{0} \right) = f\frac{\left( R_{0} - \text{rcosω} \right)}{r_{0}^{3}}, где f=6.67×10-11м3⋅кг-1⋅с-2f = 6.67 × 10^{-11} м^3⋅кг^{-1}⋅с^{-2} — гравитационная постоянная, R0,φ0,λ0R_0, φ_0, λ_0 — координаты точки измерений, r,φ,λr, φ, λ — координаты источника, r0=R02+r2-2R0rcosω r_{0} = \sqrt{R_{0}^{2} + r^{2} - 2R_{0}\text{rcosω}}, ωω — угол при центре O земного шара между точкой измерений и источником: сosω=cosφ0cosφ+sinφ0sinφcos(λ0-λ)сosω=cosφ_0 cosφ+sinφ_0 sinφcos (λ_0-λ). По сути осуществляется решение линейной обратной задачи гравиразведки [11], в котором искомые значения масс (плотностей) m не используются для получения новой геологической информации, а являются лишь средством для восстановления аномалий силы тяжести u в произвольно выбранных точках пространства, а также обеспечивают высокоточный расчет трансформант v.На первом этапе аппроксимации выполняется определение значений «глубинных» точечных масс, обеспечивающих минимум среднеквадратического расхождения исходного и модельного полей во всех имеющихся точках поля (совпадающие точки областей D1 и D2 исключаются). Для этой цели формируется нормальная СЛАУ GTGα=GTu, решение которой осуществляется методом Холецкого [12]. Эта система всегда совместна, матрица ее коэффициентов GTG обладает свойством симметричности и положительной определенности (т. н. матрица SPD — Simmetric and Positively Defined). Однако обусловленность ее несколько ухудшается: 𝑐𝑜𝑛𝑑𝐺𝑇𝐺=𝑐𝑜𝑛𝑑𝐺\sqrt{𝑐 𝑜 𝑛 𝑑 𝐺^ 𝑇 𝐺}= 𝑐 𝑜 𝑛 𝑑 𝐺 . Приблизительные потери в точности решения при выполнении условия h ≅ ΔL , где ΔL  — расстояние между точками задания поля, могут достигать 5–6 значащих цифр. Это существенно ниже критических значений. В частности, для компьютера автора статьи Asus Aspire VN7-791 с процессором Intel(R) Core(TM) i7-4710HQ экспериментально определенное значение машинного нуля 𝜀0≅1×10-19𝜀_ 0 ≅ 1 × 10^{-19}, т. е. при таких потерях результаты ориентировочно будут содержать не менее 12 значащих цифр.На втором этапе гравитационный эффект разности исходных значений поля и поля «глубинных» источников в пределах области D1 аппроксимируется системой «приповерхностных» точечных масс, значения которых определяются путем решения СЛАУ итерационным методом наискорейшего градиентного спуска (НГС) [12]. В процессе решения используется предварительно рассчитанный массив коэффициентов системы G ̃, хранящийся в оперативной памяти (RAM) компьютера. Он представляет собой разреженную матрицу, сформированную с учетом резкого уменьшения амплитуды гравитационного поля по мере удаления от точечного источника, которая на расстояниях в (15–25) h принимается равной нулю. Для гридированных данных о поле {uklu_{kl}}, k=1,m¯,l=1,n¯k = \overset{\overline{}}{1,m}, l = \overset{\overline{}}{1,n}  определение косинуса угла при центре Земли между точкой измерений поля и точечной массой сosω=cosφ0cosφ+sinφ0sinφcos(λ0-λ)сosω = cosφ_0cosφ + sinφ_0sinφ cos(λ_0 - λ)  можно существенно ускорить. Для этого достаточно однократно рассчитать mm значений cosφcosφ и nn значений cos(λ0-λ)cos(λ_0 - λ), а затем сохранить их в одномерных массивах вещественных чисел [13].Выбор относительных глубин h источников, равных шагу сети задания поля по меридиану для двух цифровых моделей поля VR, обеспечивает достаточно хорошую обусловленность СЛАУ. Это влечет за собой высокую скорость сходимости для метода НГС. По оценке авторов доклада данный принцип погружения источников не может использоваться лишь для площадей, расположенных в полярных областях Земли, за пределами широты 70°. После определения значений всех точечных масс для вычисления трансформант v в формуле (1) могут применяться различные линейные операторы T. Пересчет гравитационного поля VR на фрагмент внешней по отношению к Земле сферической поверхности с радиусом R' =  R0 + H, где H — нормальная высота, в данном случае отвечает приведению поля VZ на горизонтальную плоскость z = const для «плоской Земли». Данный алгоритм реализован в программе TRANSF_VR [14], которая позволяет осуществлять вычисление первой радиальной производной гравитационного потенциала VR, гравитационного потенциала V, его вторых производных VRλV_{Rλ}, VRϕV_{Rϕ} и модуля горизонтального градиента Gr на земной поверхности и на заданной высоте H (рис. 4). Для корректного использования производных 𝑉𝑅𝜆,𝑉𝑅𝜙𝑉_{𝑅𝜆}, 𝑉_{𝑅𝜙}  при вычислении горизонтального градиента аномалий силы тяжести Gr необходимо дополнительно использовать коэффициенты Ламе для сферической системы координат τ1=1⁄(R0sin⁡φ)τ_1=1⁄(R_{0} sin⁡φ) и τ2=1⁄R0τ_2=1⁄R_0 : Gr=(τ1VRλ2+τ2VRφ2)0.5Gr=(τ_1 V_{Rλ}^2+τ_2 V_{Rφ}^2 )^0.5.Рисунок 4. Главное окно программы TRANSF_VRТрансформация аномалий силы тяжести Сибирской платформы и обрамляющих ее структурВ качестве примера работы программы TRANSF_VR возьмем данные глобальной модели гравитационного поля EIGEN-GRGS.RL04.MEAN-FIELD, полученной на основе данных спутниковых миссий GRACE и SLR в 2019 г. Используется цифровая модель аномалий силы тяжести в полной редукции Буге с плотностью промежуточного слоя 2.67 г/см3 для территории, ограниченной 52°–72° с. ш., 84°–132° в. д.Область D1 размером ~~ 5.15 млн. км2 охватывает всю Сибирскую платформу и структуры ее обрамления. Шаг сети — 0.1°, размерность GRID-модели — 201 × 486 (N = 97 686 точек). Источником информации о рельефе земной поверхности являлась модель ETOPO1. Диапазон изменения аномалий силы тяжести ~~ 387 мГал, высотные отметки H лежат в пределах от - 19 до 2868 м. Эквивалентные источники располагаются на 11.13 км ниже земной поверхности. Поле в области D2 с координатами 5°–80° с. ш., 15°–180° в. д. задано с шагом 5°. Относительная глубина h точечных масс составляет 550 км (рис. 5).После выполнения 50 итераций методом НГС была достигнута точность аппроксимации аномалий силы тяжести ~~ 0.17 мГал в евклидовой метрике F2, что превосходит точность исходных материалов [15]. Время решения задачи составило 758 с., вычисления проводились на компьютере среднего класса OptiPlex 7070, при этом около 50 % затрат машинного времени занял ввод/вывод исходных данных и результатов расчета.Рисунок 5. Совмещенные карты изоаномал гравитационного поля для областей  D1 и D2: красный контур — территория исследований D1 После определения параметров двух разноглубинных систем точечных масс для вычисления трансформант v могут применяться различные линейные операторы T. Пересчет гравитационного поля VR  на фрагмент внешней по отношению к Земле сферической поверхности с радиусом R' = R0 + H, где H — нормальная высота, в данном случае отвечает приведению поля VZ на горизонтальную плоскость z = const для «плоской Земли». Аналогом вертикальной производной силы тяжести VZZ является вторая радиальная производная силы тяжести VRR. Соответствующую аналогию можно также провести между вторыми производными VZX, VZY и VRλ, VRφ. Вторые производные гравитационного потенциала характеризуют скорость изменения силы тяжести в пространстве и отражают локальные особенности геологического строения. Для корректного использования производных VRλ , VRφ при вычислении горизонтального градиента аномалий силы тяжести Gr необходимо дополнительно использовать коэффициенты Ламе для сферической системы координат 𝜏1=1⁄(𝑅0𝑠𝑖𝑛⁡𝜑) 𝜏_ 1 = 1 ⁄ ( 𝑅_ 0 𝑠 𝑖 𝑛 ⁡ 𝜑 ) 𝜏1 = 1 ⁄ ( 𝑅_ 0 𝑠 𝑖 𝑛 ⁡ 𝜑 ) и 𝜏2=1⁄𝑅0𝜏_ 2 = 1 ⁄ 𝑅_ 0 : 𝐺𝑟=(𝜏1𝑉𝑅𝜆2+𝜏2𝑉𝑅𝜑2)0.5𝐺 𝑟 = ( 𝜏_ 1 𝑉_ {𝑅 𝜆}^ 2 + 𝜏_ 2 𝑉_ {𝑅 𝜑}^ 2 )^ 0.5.Сопоставим представленные на 50-й сессии семинара им. Д. Г. Успенского — В. Н. Страхова результаты пересчета гравитационного поля VR в верхнее полупространство на высоту Н = 250 км [16] c новыми данными, полученными при двух уровнях размещения эквивалентных источников. Для графического отображения данных (рис. 6), как и ранее, использована равноугольная нормальная коническая проекция (НКП), в которой меридианы изображаются прямыми, сходящимися в общей точке под углами, пропорциональными разности долгот, а параллели изображаются концентрическими окружностями, имеющими общий центр в точке пересечения меридианов [17]. Параметры НКП: главная параллель 62°, центральный меридиан 108°. Как можно заметить, недостаток суммарной массы источников на периферии территории привел к появлению артефактов — замыканию изолиний вблизи ее южной границы (рис. 6а). Наличие «глубинных» масс повлекло за собой существенное изменение амплитудных характеристик аномалий силы тяжести при пересчете на высоту (рис. 6б).В системе КОСКАД 3D [18] выполнялся интеллектуальный анализ результативного набора трансформант с использованием вероятностно-статистических методов. С целью трассирования аномалий гравитационного поля на высоте Н = 25 км различных энергий и различного направления применялась оригинальная модификация одномерной адаптивной фильтрации. Выделенные оси положительных и отрицательных аномалий, связанные преимущественно с дизъюнктивной тектоникой, показаны на рисунке 7а. С помощью двухмерной энергетической фильтрации были выделены различные составляющие этого поля, наиболее энергоемкая из них (тренд) приведена на рисунке 7б. Морфология этой составляющей отражает наиболее крупные тектонические элементы Сибирской платформы и неоднородное по плотности строение верхней мантии под ней.Рисунок 6. Карты изоаномал силы тяжести VR в редукции Буге на высоте 250 км, полученные путем пересчета поля в верхнее полупространство при одном (а) и двух (б) уровнях размещения точечных массРисунок 7. Результаты автоматического трассирования осей аномалий (а) и наиболее энергоемкая составляющая гравитационного поля VR (б)Объемное распределение эффективного параметра, пропорционального плотности горных пород, построено в системе КОСКАД 3D с применением статистических, спектрально-корреляционных методов и алгоритма адаптивной фильтрации в окне живой формы на основе преобразования Б. А. Андреева. Блок-диаграмма квазиплотности (рис. 8) представляет собой результат перехода от 2D- к 3D-геоизображениям на основе аналитических компьютерных процедур логико-математической генерализации. Последовательное повышение уровня генерализации обеспечивает проявление на геоизображениях элементов все более крупных геосистем, что способствует выявлению качественно новой информации и закономерностей. Построенная блок-диаграмма характеризует основные особенности глубинного строения Сибирской платформы и обрамляющих ее структур, в т. ч. слабо проявленные в наблюденном гравитационном поле.Результаты свертки информации методом компонентного анализа по 7 признакам (гравитационное поле на высоте 25 км, модуль горизонтального градиента на той же высоте, 2-я радиальная производная гравитационного потенциала на высоте 75 км, 4 компоненты разложения поля методом энергетической фильтрации) приведены на рисунке 9. Пространственное распределение 1-й главной компоненты имеет отчетливо выраженную взаимосвязь с известными тектоническими структурами Сибирской платформы — Енисей-Хатангским перикратонным прогибом, Тунгусской и Вилюйской синеклизами, Алдано-Становым щитом и др.Рисунок 8. Блок-диаграмма относительного распределения плотности в литосфереРисунок 9. Карта 1-й главной компоненты, включающей в себя 7 признаков — трансформант гравитационного поляВыводыРазработана компьютерная технология трансформации аномалий силы тяжести в пределах больших территорий, в которой используются «квазиэллипсоидальная» модель Земли — сфера Каврайского и геодезические координаты точек задания поля. При вычислениях используются две телескопированные цифровые модели гравитационного поля с разным шагом сети при двух уровнях размещения точечных масс — «глубинном» и «приповерхностном». Это позволяет уменьшить негативные эффекты, проявляющиеся при одной глубине размещения эквивалентных источников и обеспечивает высокую скорость вычислительного процесса. Результативный набор трансформант, представленных в цифровой форме, может использоваться для интеллектуального анализа в системе КОСКАД 3D методами интерпретационной томографии, автоматизированной классификации и распознавания образов. Также предварительно может осуществляться геологическое редуцирование интерпретируемого гравитационного поля.БлагодарностиИсследование выполнено при финансовой поддержке Министерства науки и высшего образования РФ в рамках государственного задания (рег. номер НИОКТР: 124020500054-3).