Физика Земли, 2023, № 2, стр. 62-77

ВВЕДЕНИЕ

Изучение строения верхней части литосферы (ВЧЛ) гравиметрическими методами и построение региональных геофизических моделей опирается на решение следующей задачи: найти распределение модельной плотности, которое удовлетворяет наблюденному гравитационному полю. Интегральная формула прямой задачи – это оператор “суммы”, заданный на множестве плотностных элементов дискретного разбиения. И найти плотности каждого отдельного элемента только по их совокупному гравитационному эффекту теоретически невозможно. В общей постановке решение обратной задачи гравиметрии является не только неоднозначным, но и неустойчивым. И чем выше размерность задачи, тем она менее устойчива и тем настоятельнее требуется применение специальных регуляризирующих алгоритмов. Регуляризация обеспечивает устойчивость решения обратной задачи, но не дает его единственности. Множество решений регуляризирующих алгоритмов образуют эквивалентное семейство по параметру регуляризации. Поэтому, вопрос о выборе единственного решения нельзя рассматривать только с позиций общей математической теории; необходимо привлекать “сторонние” плотностные модели начального приближения. Начальная модель, построенная посредством послойной интерполяции сеточных матриц с двумерных плотностных разрезов, сохраняет основные структурные элементы глубинного строения в окрестности сейсмических профилей и транслирует их подобие в межпрофильное пространство [Мартышко и др., 2012; 2016а]. Разделение плотностного (сеточного) параллелепипеда начальной модели на семейство горизонтальных слоев реализует алгоритм поиска слабоединственных решений обратной задачи гравиметрии для трехмерной избыточной плотности на множестве корректности двумерных мультипликативных добавок [Мартышко и др., 2013; 2016б]. По каждому горизонтальному слою вычисляются малые (латеральные) отклонения плотности от ее начального распределения. Это дает возможность из семейства послойных эквивалентов избыточных плотностей выбрать единственное частное решение с минимальной нормой.

Настоящая работа является продолжением исследований по созданию устойчивых математических методов трехмерного гравитационного моделирования применительно к решению практических задач и построению трехмерной геофизической модели глубинного строения ВЧЛ для территории Северного и Среднего Урала. Новые систематизированные данные позволяют расширить (в направлении на юг) пространственную область количественной интерпретации гравитационных полей с построением сеточной цифровой модели трехмерного распределения плотности в земной коре и верхней мантии. Настоящие результаты являются обобщением ранее опубликованных исследований по 4-м сейсмическим профилям [Мартышко и др., 2022] и завершают построение пространственной плотностной модели на основе 10 плотностных разрезов и устойчивой методики решения обратной задачи гравиметрии в классе слабо единственных эквивалентов.

СИСТЕМАТИЗИРОВАННЫЕ ДАННЫЕ ДЛЯ АЛГОРИТМА ТРЕХМЕРНЫХ ВЫЧИСЛЕНИЙ

Обратная задача вычисления трехмерного распределения плотности по аномалиям гравитационного поля опирается на объемную плотностную модель начального приближения, реконструированную вдоль сейсмических профилей. Современные представления о развитии Урала и его платформенного обрамления учитывались в рамках существующих структурных схем тектонического районирования отдельных геологических провинций. Для Уральской складчатой системы (УСС) наиболее содержательной является схематическая карта Урала М: 1 : 2 500 000 под редакцией И.Д. Соболева [Соболев, 1968].

На рис. 1 приведен фрагмент гравитационного поля в полной редукции Буге, взятый с комбинированной глобальной модели гравитационного поля XGM2019e_2159_GA [Zingerle et al., 2020] и сконвертированный в координатах картографической проекции Гаусса–Крюгера. Положение региональных сейсмических профилей привязано к фрагменту карты поля; структурные схемы тектонического районирования использованы для верификации результатов количественной интерпретации гравиметрических данных.

ДВУМЕРНЫЕ СКОРОСТНЫЕ РАЗРЕЗЫ

Информационная база для построения плотностной модели исследуемой территории включает три составляющие: поля времен по профилям ГСЗ и соответствующие им двумерные скоростные разрезы земной коры, эмпирическую (предварительную) корреляционную зависимость “скорость–плотность” и цифровые карты площадных аномалий гравитационного поля. Для модели ВЧЛ миллионного масштаба данные должны быть заданы на сетке с шагом не более 1 км. Скоростные разрезы, построенные по профилям ГСЗ, и их двумерные плотностные аналоги не смогут воспроизвести детали пространственных объектов только по видимым проективным сечениям [Егоркин, 1991; Атлас…, 2013]; необходимую трехмерную информацию можно извлечь из карт гравитационных и магнитных аномалий соответствующего масштаба [Павленкова, Романюк, 1991].

Для построения трехмерной сейсмоплотностной модели Среднеуральского сегмента в пределах заявленной территории (56°–60° с.ш. и 54°– 60° в.д.) были использованы данные сейсмических исследований, охватывающих существенно большую площадь – градусную трапецию (53–65° с.ш., 48–72° в.д.). Такой формат задания исходных данных предполагает построение пространственного сейсмического каркаса оцифрованных скоростных разрезов для расширенной территории с последующим замыканием контура интерполяции плотностной модели на меньшую площадь (рис. 1).

На территории Среднеуральского сегмента располагается десять региональных сейсмических профилей, восемь из которых выполнены Баженовской геофизической экспедицией ПГО “Уралгеология” (три профиля – при участии Института геофизики УрО РАН), и двух геотраверсов Центра региональных геофизических и геоэкологических исследований имени В.В. Федынского (Центра “ГЕОН”). По фондовым материалам для всех сейсмических профилей выполнена обработка материалов ГСЗ по методу двумерной сейсмической томографии [Мишенькина и др., 1983]. Для трех профилей Баженовской геофизической экспедиции (Тараташский, Красноленинский и Северная Сосьва-Ялуторовск) были переинтерпретированы ранее построенные двумерные поля времен первых вступлений продольных Р-волн. Для остальных профилей: Вижай–Нижняя Тура–Орск, Свердловский, Гранит -Рубин-2, а также для объединения двух профилей (Красноуральский + Ханты-Мансийский) и фрагмента геотраверса Рубин-1 поля времен были построены заново. Эти поля времен для десяти разрезов ГСЗ легли в основу сеточных (градиентных) двумерных скоростных моделей, рассчитанных на всю мощность земной коры по авторской программе алгоритма “Inverse”. Дополнительно, по фрагментам двух геотраверсов центра “ГЕОН” (Рубин-1 и Рубин-2) слоисто-блоковые модели были переформатированы в двумерные матрицы сеточных функций. Таким образом были перестроены скоростные разрезы для западного участка геотраверса Рубин-1 (от 48° до 55° в.д.) и восточного участка Рубин-2 (от 67° до 72° в.д.). Все разрезы дополнены границей Мохо (М), положение которой определено по скоростным уровням (7.75–8.25) км/с, а затем откорректировано по результатам интерпретации имеющихся данных по отраженным, обменным и головным волнам. Погрешность построения скоростных разрезов оценивалась при сопоставлении их в местах пересечения (или сближения) профилей. Общая протяженность профилей ГСЗ более 7500 км.

СЕТОЧНЫЕ ПЛОТНОСТНЫЕ МОДЕЛИ ПО СКОРОСТНЫМ РАЗРЕЗАМ

Прагматическая цель построения плотностных моделей – проверка соответствия скоростных параметров аномальному гравитационному полю по совокупности имеющихся разрезов. Преобразование скоростных параметров земной коры в плотностную матрицу и технология построения блочной структуры верхней мантии вдоль Ханты-Мансийского профиля изложена в работе [Мартышко и др., 2022]. Итоговый результат схематично изображен на рис. 2.

Обработка скоростных матриц и построение градиентной слоисто-блочной модели распределения плотности по анонсированным выше профилям выполняются по той же схеме. Перечислим ее основные этапы.

1. По полям времен продольных волн в первых вступлениях строится сеточный скоростной разрез земной коры [Крылов и др., 1993] . По кусочно-линейной регрессионной зависимости “скорость–плотность” рассчитывается массив предварительных значений плотности на сетке той же размерности [Дружинин и др., 1982; Семенов, 1993]. Модель дополняется мантийным слоем с плотностью ${{\sigma }_{M}} = {\text{3}}{\text{.3 г/с}}{{{\text{м}}}^{{\text{3}}}}$ до глубины регионального уровня изостатической компенсации 80 км. Соответствие скоростных и плотностных параметров вынесено на шкалы петрофизической колонки (см. рис. 2.1 ).

2. По аномалиям литостатического давления строится блочный разрез верхней мантии. Изостатические поправки к плотности – функция-компенсатор ${{\Sigma }}\left( x \right)$, позволяют наметить контуры мантийных блоков и оценить в них значения плотности (см. рис. 2.2 ).

3. Рассчитываются значения одномерной гидростатической плотности как средние по горизонтальным слоям плотностной модели. Вмещающая среда за пределами профиля заполняется массами с гидростатической плотностью. Относительно последней вычисляются избыточные плотности разреза. При этом гравитационный эффект от вмещающей среды равен нулю, а разностное поле самого разреза будет приведено к нулевому уровню.

4. Составляется функционал модельного гравитационного поля, в котором неизвестные элементы плотностной матрицы выражаются через интервальные значения коэффициентов кусочно-линейной зависимости “скорость–плотность” в земной коре. Формируется искомый вектор неизвестных коэффициентов и плотности мантийных блоков. По каждому скоростному интервалу уточняющая заверка распределения плотности выполняется в процессе решения 2D-линейной обратной задачи гравиметрии. Подбираются не плотностные элементы сеточных матриц, число которых порядка 4 × 10⁴, а только четыре коэффициента кусочно-линейной регрессии и плотности мантийных блоков. Этого достаточно, чтобы уточнить матрицу избыточных плотностей сводной модели “кора–мантия” (см. рис. 2.3 ).

5. Подобранные сеточные плотности по разрезам до глубины 80 км сопоставляются с сеточным массивом исходных данных по скоростям продольных волн в коре и на границе “кора–мантия”. Связь двух однотипных числовых массивов позволяет уточнить исходную петрофизическую зависимость “плотность–скорость” и, в качестве альтернативы, определить поинтервальные значения материальных коэффициентов “А” и “В” по результатам численного моделирования.

В таблице приведены уточненные коэффициенты кусочно-линейной зависимости “плотность–скорость” по каждому из профилей и мера несогласия (среднеквадратичное отклонение RMS) модельного и наблюденного полей. Оценки даны как для плотностных разрезов с однородной мантией (начальное приближение), так и для подобранных разрезов с мантией блочной структуры.

На рис. 3 представлены двумерные сейсмоплотностные модели земной коры и верхней мантии до глубины изостатической компенсации 80 км. Расчетное гравитационное поле моделей сопоставлено с наблюденными аномалиями Буге. Морфологическое сходство двух кривых вполне приемлемо на качественном уровне интерпретации.

НАЧАЛЬНАЯ МОДЕЛЬ ИНТЕРПОЛИРОВАННОЙ ПЛОТНОСТИ

Двумерные скоростные и соответствующие им плотностные разрезы по десяти сейсмическим профилям формируют исходную базу данных трехмерной плотностной модели начального приближения. Планшет профильных данных покрывает площадь, существенно большую заявленной территории – градусной трапеции (53–65° с.ш., 48–72° в.д.). Увеличение длины сейсмических профилей предполагает построение пространственного каркаса оцифрованных разрезов для расширенной территории с последующим замыканием контура интерполяции на требуемый сегмент. Это позволило снизить влияние градиентных зон из-за разрыва интерполированной плотности на границах планшета.

Двумерные сеточные матрицы с плотностных разрезов форматируются в координатах цифрового планшета карты аномалий гравитационного поля. В предложенной схеме допустимо использование разрезов различной детальности (масштаба) и глубины. Граничные файлы образующих кривых задают положение каждого цифрового фрагмента в площадном варианте. Таким образом учитывается взаимное расположение в пространстве и происходит переход от 2D-массива координат вертикальных сечений к 3D-сеточному массиву координат объемной модели, по которому и восстанавливается макет пространственного сейсмического каркаса (рис. 4).

Недостающие данные по плотностям в межпрофильном пространстве заполняются интерполированными значениями плотности с опорных разрезов трехмерного каркаса. Интерполяция выполняется по отдельным горизонтальным слоям: мощности слоев определяются выбранным нами шагом дискретизации (100 м) по оси глубин; в горизонтальных плоскостях шаг интерполяции 1 км. При этом учитывается схема тектонического районирования и, по возможности, присущая зональность гравитационного поля. Данные за пределами граничного контура не учитываются. Результатом работы программы интерполяции является четырехмерный массив послойных сеточных файлов (x, y, z, σ₀) и XML-файл, в котором каждому слою сопоставлена глубина, на которой он расположен:

Средние значения плотности по каждому слою приравнены к распределению одномерной гидростатической плотности, зависящей только от глубины:

Плотностной параллелепипед интерполированной плотности ${{\sigma }_{0}}\left( {x,y,z} \right)$ дополняется до 3D-модели неоднородного горизонтального пласта бесконечного простирания. Для этого внешнее законтурное пространство заполняется массами с одномерной гидростатической плотностью ${{\sigma }^{0}}(z)$. Относительно последней вычисляются избыточные плотности. Таким образом, параллелепипед избыточных плотностей будет находится во вмещающей среде с нулевым гравитационным эффектом [Ладовский и др., 2017].

Поле модели вычисляется по трехмерной матрице значений избыточной плотности для пласта бесконечного простирания с выделением конечного фрагмента в пределах целевой области задания наблюденного поля. В том же формате, что и наблюденное поле, вычисляется и разностное поле. Программа вычислений “3D_σ CALC” составлена на основе “быстрых” алгоритмов и обладает высокой производительностью [Мартышко др., 2013]. Так, для модели целевой области с разбиением 500 × 750 × 800 элементов время расчета поля на сетке 500 × 750 на одном графическом процессоре NVidia Titan Black с тактовой частотой до 3 ГГц (при параллельной программной реализации) не превышает 18 мин. На рис. 5 показан фрагмент общего решения трехмерной прямой задачи гравиметрии в пределах заданного участка исследуемой территории.

Выявленное морфологическое отличие модельного (трехмерного) поля от наблюденного подчеркивает несогласие исходных 2D- и построенных 3D-моделей интерполированной плотности. Амплитуда изменений разностного поля служит критерием точности построения сейсмоплотностной модели начального приближения и сопоставима с амплитудой наблюденного поля. Вариант подобранной плотностной модели надежнее строить по разработанной нами технологии трехмерного сеточного гравитационного моделирования. В разностное поле уже включена информация по профилям о плотностной модели начального приближения. По этой причине именно числовой массив разностного поля принимается за исходный массив входных данных при решении линейной обратной задачи гравиметрии на завершающем этапе построения трехмерных плотностных моделей по аномалиям наблюденного гравитационного поля.

РАЗДЕЛЕНИЕ РАЗНОСТНОГО ПОЛЯ ПО ГОРИЗОНТАЛЬНЫМ СЛОЯМ ИСТОЧНИКОВ

Обратная задача гравиметрии – вычисление трехмерной плотности неоднородной области по заданным на множестве внешних точек значениям гравитационного поля, является классическим примером некорректной задачи: в общей постановке ее решение неединственно и неустойчиво зависит от исходных данных. Поэтому необходимо искать решения на практически содержательных множествах корректности (классах) модельных распределений плотности. Так, из класса слабоединственных эквивалентов можно выбрать устойчивое решение обратной задачи гравиметрии при условии, что размерность неизвестной плотности совпадает с размерностью внешнего поля (зависит от одних и тех же переменных) [Новоселицкий, 1965]. С точки зрения дискретных матричных моделей наиболее приемлемым оказался вариант восстановления латерально изменчивой плотности в горизонтальном слое. Но для этого необходимо из совокупности наблюденных данных выделить поля отдельных слоев. Метод сглаживания наблюденного поля при пересчете на различные высоты и последующего продолжения “вниз” может использоваться для фильтрации и послойного разделения аномалий от источников, локализованных в горизонтальных слоях на соответствующих глубинах [Martyshko et al., 2021].

Результат фильтрации разностного поля с использованием методов аналитического продолжения на глубину показан на рис. 6. Слева приведен пример разделения поля на составляющие от четырех горизонтальных слоев земной коры с указанием соответствующих границ по глубине и значениями амплитуд разделенных аномалий. Справа дана “рабочая” схема в виде параллелепипеда разделенных полей по разноглубинным срезам с шагом 1 км, которая далее будет использована, как база входных данных при решении обратной задачи гравиметрии для послойного восстановления трехмерной избыточной плотности.

Результат скомпилирован из работы авторов [Martyshko et al., 2021], где подробно излагается вычислительный алгоритм разделения полей при их пересчете на несколько высот и последующем аналитическом продолжении на глубину с выбором параметра адаптивной регуляризации.

ЛИНЕЙНАЯ ОБРАТНАЯ ЗАДАЧА ГРАВИМЕТРИИ ДЛЯ МОДЕЛИ РАЗДЕЛЕННЫХ СЛОЕВ

Вычисление трехмерной плотности $\sigma \left( {x,y,z} \right)$ неоднородной области по заданным на множестве внешних точек значениям поля $g\left( {x,y} \right)$ (обратная задача гравиметрии) реализуется на основе решения операторного уравнения первого рода. Для горизонтальных слоев с двумерным распределением плотности (параметрическом классе корректности двумерных распределений) устойчивое решение обратной задачи гравиметрии существует и однозначно соотносится с пространственными аномалиями исходного (наблюденного) гравитационного поля. Алгоритм поиска устойчивых решений опирается на принцип локализации – варианте локально-одномерного проектирования трехмерного поля на горизонтальную плоскость с последующей минимизацией нормы невязки для плоско параллельного слоя расчетного и наблюденного полей.

Зададим геометрию исходного плотностного параллелепипеда ${{D}_{H}}:$

и элементов горизонтальных пластов ${{D}_{K}}$ его послойного разделения

Разностное поле $\Delta g = g\left( {x,y} \right) - {{g}_{0}}\left( {x,y} \right)$ вычисляется как отклонение наблюденного поля от расчетного поля модели начального приближения. Эта разность позиционируется как целевая функция подбора избыточной плотности $\Delta \sigma = \sigma \left( {x,y,z} \right) - {{\sigma }_{0}}\left( {x,y,z} \right)$.

Устойчивость решения обратной задачи гравиметрии базируется на представлении избыточной трехмерной плотности $\Delta \sigma $ в виде мультипликативной функции (произведения одномерной гидростатической плотности начальной модели на переменную по латерали корректирующую добавку в каждом слое) [Мартышко и др., 2013; Martyshko et al., 2018]:

где – латеральная корректирующая добавка избыточной плотности в k-ом слое; в выбранном классе двумерных плотностей решение ${{\Phi }}\left( {x,y} \right)$ единственно.

Модель начального приближения разбивалась на несколько горизонтальных слоев ${{D}_{k}}$ по интервалам постоянства гидростатической плотности ${{\sigma }^{0}}\left( z \right)$ (здесь ${{z}_{k}} - {{z}_{{k - 1}}} = 1\,\,{\text{км}}$). Каждому слою ставились в соответствие разделенные аномалии ${{\delta }}{{g}_{k}}$ повысотных трансформант разностного поля $\Delta g = \sum\nolimits_k {\delta {{g}_{k}}} $ (см. рис. 6); избыточные плотности $\delta {{\sigma }_{k}} = {{\sigma }^{0}}\left( z \right){{\Phi }_{k}}\left( {x,y} \right)$ в выделенных слоях однозначно вычисляются по разделенным полям $\delta {{g}_{k}}$ (принцип локализации латеральных вариаций). Функция ${{{{\Phi }}}_{к}}\left( {x,y} \right)$ в k-ом слое удовлетворяет интегральному уравнению Фредгольма 1-го рода. Поле начальной модели $U\left( {\xi ,\eta ,\zeta } \right) = \sum\nolimits_k {{{U}_{k}}} $ вычисляется по каждому отдельному слою ${{D}_{K}}$ (${{\gamma }}$ – гравитационная постоянная):

Эта функция задает начальную разность итеративного подбора при решении интегрального уравнения:

Таким образом, задача вычисления корректирующей добавки к распределению плотности для одного слоя плотности сводится к двумерному случаю и может решаться независимо от других слоев. Это обеспечивает единственность решения уравнения для слоя ${{{{\Phi }}}_{k}}\left( {x,y} \right)$ и, как следствие, для всего параллелепипеда ${{\Phi }}\left( {x,y} \right) = {{ \cup }_{{\left( k \right)}}}{{{{\Phi }}}_{k}}$ на основе устойчивого численного алгоритма инверсии.

Решение обратной задачи гравиметрии в формате сеточных функций находится на основе принципа локализации. В работе [Мартышко и др., 2016б; Martyshko et al., 2021] для уравнения Фредгольма предложен соответствующий итерационный алгоритм для нахождения значений ${{\Phi }_{k}}\left( {x,y} \right)$, минимизирующий невязку $\delta g$ наблюденного (разностного) поля $\delta {{g}_{K}}\left( {\xi ,\eta ,\zeta } \right)$ и модельного ${{U}_{K}}\left( {\xi ,\eta ,\zeta } \right)$ полей по каждому отдельному слою и всей модели в целом. Сопоставление нормы невязки поля для ряда последовательных приближений показывает, что уже после 20-й итерации амплитуда невыбранных остатков разностного поля уменьшается по сравнению с исходным более чем в 40 раз. Для исключения влияния краевых эффектов все расчеты выполнялись для расширенной территории. Инверсия поля реализована для территории (53–65° с.ш., 48–72° в.д) на сетке (ШИРОТА × ДОЛГОТА × ГЛУБИНА = = 1236 × 1314 × 80 км); размер ячейки 1 × 1 × 1 км. В этом случае одна итерация для 80-и слоев считается 2 мин. Избыточная плотность послойно разделенной модели восстанавливается мультипликативной функцией гидростатической плотности ${{\sigma }^{0}}\left( z \right),~~z \in \left( {{{z}_{{k - 1}}},~~{{z}_{k}}~} \right)$ и латеральной корректирующей добавки ${{\Phi }_{k}}\left( {x,y} \right)$ с нулевым средним значением. Сеточные матрицы искомых решений скомпонованы в единый пространственный макет числового параллелепипеда; на его основе строится объёмная модель послойного распределения избыточной плотности. Физическая плотность модели (в абсолютной мере) является суммой подобранного (с нулевым средним) и начального распределений; в таком случае, характер изменения гидростатической плотности ${{\sigma }^{0}}\left( z \right)$ по глубине остается практически неизменным:

На рис. 7 построено искомое распределение послойно-подобранной избыточной плотности для сеточного параллелепипеда 500 × 750 × 80 км, что соответствует размерам целевой области. Финальная модель изображена в пределах градусной трапеции 56°–60° с.ш. и 54°–66° в.д. (см. рис. 1а).

Сеточные решения обратной задачи гравиметрии, сконвертированные в формате числового параллелепипеда, дают наглядное представление о распределении плотности в пространстве объемной модели. Из общего трехмерного массива подобранных плотностей $\sigma \left( {x,y,z} \right)$ возможно извлечение любых массивов меньшей размерности для трассировки вертикальных сечений, построения карт горизонтальных срезов или структурных карт сейсмо-плотностных границ (с привязкой к схемам тектонического районирования). Последнее весьма удобно для всестороннего анализа и последующей геолого-геофизической интерпретации результатов математического моделирования [Мартышко и др., 2016а]. Так например, на рис. 8 даны два варианта плотностных сечений вдоль объединенного профиля “Гранит–Рубин 2” для начальной и результирующей плотностной модели, извлеченные из послойного набора сеточных матриц по направлению образующих профильных кривых.

Как видно, при переходе от начальной модели интерполированной плотности к результирующей модели, структурные особенности разрезов, в целом, сохранены, но общая интегральная масса коры, а следовательно, и абсолютный уровень фонового поля существенно меняются. Подчеркнем две интересные (и геологическизначимые) детали. В разрезе приповерхностных слоев мульдообразной формы выделяется приповерхностный тяжелый блок в районе 60° (район Уральской сверхглубокой скважины СГ-4); возрастает (сравнительно с сейсмоплотностным разрезом) мощность разуплотненных пород осадочного палеозоя и мезо кайнозоя на Восточном склоне Урала и Западной Сибири. Вычисление поля для всех разрезов подобранной модели и сравнение их с сейсмоплотностными разрезами начальной модели указали бы, в дальнейшем, на зоны некондиционного материала сейсмических зондирований.

Практическим и весьма наглядным итогом количественной интерпретации гравитационных аномалий является сравнительный анализ одномерных характеристик, извлеченных из трехмерных массивов исходных и рассчитанных плотностей. На рис. 9 приведены зависимости от глубины минимальных, максимальных и средних значений плотности начальной и подобранной трехмерных моделей. Для каждой глубины они вычислялись по множеству значений плотности в точках сечения модели горизонтальной плоскостью с аппликатой, равной это глубине.

Сравнивая эти одномерные распределения, нетрудно заметить, что поправки в значения плотности за гравитационное моделирование не достигают глубинных горизонтов в низах коры и верхней мантии. Проявляется аномальный гравитационный эффект от источников, сосредоточенных в верхней (порядка 30 км) толще. И здесь мы сталкиваемся с неким аналогом “гравиактивного слоя”, содержательным распределением плотности в котором можно удовлетворить аномальному полю.

ЗАКЛЮЧЕНИЕ

Для корректного решения обратной задачи гравиметрии и построения объемных плотностных моделей литосферы по сейсмическим профилям потребовалась предварительная подготовка и дополнение исходных 2D-данных, их специальная конвертацая в 3D-формат послойных сеточных функций и создание пространственного макета адекватной математической модели. Методическое сопровождение алгоритма трехмерного плотностного моделирования включает следующие этапы.

1. Томографический скоростной разрез (двумерная матрица) трансформируется в свой плотностной аналог по коэффициентам кусочно-линейной петрофизической зависимости “скорость–плотность”, изначально полученной из эксперимента. Плотностной разрез замыкается снизу горизонтальной границей изостатической компенсации (по Дж. Пратту) для неоднородной мантии. По аномалиям литостатического давления намечаются контуры мантийных блоков.

2. Для плотностного разреза решается обратная линейная задача гравиметрии, в результате чего уточняются плотности мантийных блоков и значения коэффициентов кусочно-линейной функции “скорость–плотность”. Для устранения краевых эффектов законтурное пространство за пределами разреза заполняется внешними массами с одномерной гидростатической плотностью.

3. Двумерные плотностные разрезы вдоль десяти профилей ГСЗ Северо – и Средне Уральского сегмента пересчитываются в единую систему пространственных координат в пределах градусной трапеции задания исходных данных. В этих координатах строится пространственный каркас трехмерной плотностной модели. Недостающие данные по плотностям в межпрофильном пространстве заполняются интерполированными значениями плотности с опорных разрезов трехмерного каркаса. Послойное заполнение сеточных матриц числового параллелепипеда сохраняет подобранные 2D-плотностные параметры в сечении сейсмических профилей.

4. Пространственная модель послойно-интерполированной плотности принимается за модель начального приближения. Относительно нее рассчитывается избыточная плотность и разностное поле целевой функции уточняющего подбора. Для слоистого параллелепипеда избыточной плотности гравитационный эффект законтурных (внешних) масс равен нулю. Сеточный формат редуцированной математической модели сделал возможным использование итеративной схемы послойного подбора плотностных разностей, образующих пространство корректных решений трехмерной обратной задачи гравиметрии.

5. По формулам пересчета “вверх–вниз” разностное поле разделяется на морфо подобные разнопорядковые составляющие, каждая из которых ассоциируется с модельным полем слоев редуцированной плотностной модели. Избыточные плотности в каждом слое заданы в виде мультипликативной функции (произведения гидростатической плотности на переменную по латерали корректирующую добавку). Устойчивый алгоритм послойной линейной инверсии реализован в классе двумерных корректирующих добавок с нулевым средним значением. Итеративный алгоритм последовательных вычислений избыточных плотностей применен для каждого отдельного слоя, независимо от других. Трехмерные мультипликативные плотности по слоям однозначно восстанавливаются по отфильтрованным разностям повысотных трансформант. Таким образом удалось сократить до минимума число возможных вариантов подбора разделенных полей и унаследовать геологическую содержательность плотностной модели начального приближения.

Сеточный формат задания исходных данных, “быстрые” алгоритмы решения прямой задачи гравиметрии, построение однотипных функционалов на послойных сеточных матрицах и реализация параллельных многопоточных вычислений на графических процессорах образуют основу методов количественной интерпретации геофизических полей. Все эти алгоритмические разработки предназначены для построения крупномасштабных математических моделей, как будущей основы реальных геологических проектов.