Физика Земли, 2022, № 3, стр. 12-25

1. ВВЕДЕНИЕ

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

Геолого-геофизические модели глубинного строения земной коры и верхней мантии (верхней части литосферы – ВЧЛ) являются итоговым результатом комплексной интерпретации всех имеющихся данных. При этом одним из основных индикаторов корректной реализации методики комплексирования геофизических полей служит плотность [Красовский, 1989]. Именно в ней наиболее полно отражаются петрофизические особенности неоднородного строения и литологического состава горных пород. Способ и технология расчета трехмерного распределения плотности с привязкой к двумерным данным по опорным сейсмическим разрезам закладываются в разработанную нами методику количественной интерпретации потенциальных полей с построением объемных геофизических моделей.

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

Технология построения крупномасштабных трехмерных плотностных моделей должна быть оптимальной по алгоритму вычислений, что предполагает адекватное представление исходного сейсмического материала. Авторы применили сеточный формат для построения “градиентных” скоростных и плотностных разрезов земной коры. Связь между скоростными и плотностными параметрами сеточных моделей изначально выполняется по эмпирической регрессионной зависимости “плотность–скорость”, заданной в двух скоростных интервалах для эффузивно-осадочных и кристаллических пород различного вещественного состава [Дружинин и др., 1982]. Уточняющая заверка коэффициентов искомых кусочно-линейных функций выполняется в процессе решения 2D-обратной задачи гравиметрии по профильным данным. Подобный подход применялся в работах [Страхов, Романюк, 1984; Романюк, 1995].

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

На рис. 1 приведена схема региональных профилей, совмещенная с фрагментом карты аномалий гравитационного поля с листов карты ОP 40–ОP 42 картографической проекции Гаусса–Крюгера [Zingerle et al., 2020]. Сводная схема тектонического районирования составлена на основе тектонической карты Урала [Соболев, 1969]. Пунктиром нанесены границы Зауральских структур, а именно: западная граница Ханты-Мансийского срединного поднятия [Мегакомплексы…, 1986] и восточная граница Тюменско-Кустанайского прогиба, выделенная нами по профилям “Ханты-Мансийский” и “Северная Сосьва–Ялуторовск”.

Рис. 1.

Профили ГСЗ на Северном Урале (цифры в кружочках): 1 – Красноленинский, 2 – Красноуральский и Ханты-Мансийский, 3 – Вижай–Нижняя Тура–Орск, 4 – Северная Сосьва–Ялуторовск. Структуры первого порядка по И. Д. Соболеву [Соболев, 1969]: I – Восточно-Европейская платформа, II – Предуральский прогиб, III – Западно-Уральская зона линейной складчатости, IV – Центрально-Уральское поднятие, V – Тагильско-Магнитогорский прогиб, VI – Восточно-Уральское поднятие, VII – Восточно-Уральский прогиб, VIII – Зауральское поднятие, IX – Тюменско-Кустанайский прогиб, X – Ханты-Мансийское срединное поднятие.

2. СЕТОЧНЫЕ МОДЕЛИ СКОРОСТНЫХ РАЗРЕЗОВ ГСЗ. МЕТОД СЕЙСМИЧЕСКОЙ ТОМОГРАФИИ

Результаты глубинных сейсмических исследований на профилях Северного Урала (см. рис. 1) приведены в ряде обзорных работ [Дружинин, 1983; Дружинин и др., 1990] и обобщающей монографии [Дружинин и др., 2014]. Слоисто-блоковые сейсмические разрезы земной коры построены с использованием традиционных схем интерпретации преломленных, отраженных и обменных волн. На разрезах, как правило, показаны опорные сейсмические границы, соответствующие поверхностям фундамента и мантии, а также внутрикоровые поверхности с оценкой граничных, пластовых и средних скоростей. К сожалению, значения пластовых скоростей оценены не равномерно по разрезу, а зависят от плотности систем наблюдений и от возможности непрерывной идентификации (прослеживания) головных и отраженных волн. Непрерывное поле скоростей указано лишь по верхней части некоторых разрезов, обычно, до глубины 5–10 км (на Красноленинском профиле – до 30 км). Различный тип и форма представления сейсмического материала затрудняют применение автоматизированных алгоритмов численных схем расчета для решения задачи гравиметрии. Предварительно необходимо выполнить конвертацию входных данных к единому цифровому стандарту. Наиболее приемлемым оказывается сеточный формат непрерывных скоростных разрезов, например, градиентный разрез в изолиниях скорости.

При плотностном моделировании из всей доступной сейсмической информации предпочтение отдается цифровым сеточным данным, полученным при специальной обработке профильных материалов ГСЗ. Непрерывное заполнение сеточных скоростных матриц на всю мощность земной коры достигается преобразованием дискретной системы годографов в непрерывное поле времен. По отдельным годографам дифференциальных сейсмических зондирований с помощью интерполяции строится двумерное поле времен первых вступлений, что позволяет в каждой точке профиля иметь полный набор времен пробега упругих волн на разных расстояниях (базах) от источника и значения расчетных скоростей этих волн на разных глубинах [Пузырев и др., 1975]. Послойное сканирование разрезов земной коры с построением “непрерывных” скоростных матриц получило название “метод двумерной сейсмической томографии”. Один из первых вариантов томографического метода на временных задержках, основанный на линеаризованном алгоритме решения (“Invers”) обратной кинематической задачи, был разработан под руководством С.В. Крылова в ИГ СО РАН [Мишенькина и др., 1983; Крылов и др., 1993]. Наша версия программы “Invers” выдает результат в виде сеточных скоростных матриц специального вида, предназначенных для построения континуальных моделей градиентных скоростных разрезов.

В целях равномерного представления данных о распределении скоростных параметров на всю мощность земной коры авторами выполнена переобработка уральских материалов ГСЗ по схеме алгоритма “Invers”. По первичным архивным данным построены специальные двумерные поля времен первых вступлений продольных Р-волн для профилей Баженовской геофизической экспедиции: Вижай–Нижняя Тура–Орск, Красноуральский и Ханты-Мансийский, а также переинтерпретированы ранее построенные поля по профилям Красноленинский и Северная Сосьва–Ялуторовск. Поля времен пяти разрезов ГСЗ легли в основу сеточных градиентных скоростных моделей земной коры Северного Урала в пределах градусной трапеции 58°–64° с.ш., 54°–72° в.д. (рис. 1).

Градиентные модели скоростных разрезов земной коры V_P(x, z) построены в формате сеточных функций. Верхняя граница разрезов z = 0 (км) соответствует уровню земной поверхности; нижняя – границе Мохоровичича M(x, z). Положение границы М определено по скоростным уровням (7.70–8.2) км/с и затем откорректировано по результатам интерпретации имеющихся данных по отраженным, обменным и головным волнам. На рис. 2 приведен скоростной разрез земной коры в изолиниях V_P(x, z) по объединенному профилю Красноуральский + Ханты-Мансийский. Длина годографов на профиле ограничена базами 260–280 км, соответственно, глубина проникания сейсмического луча составляет 40–50 км. Поэтому двумерная модель распределения лучевых скоростей построена только в пределах земной коры. Скоростные (и плотностные) неоднородности верхов мантии на первом этапе не входят в нашу схему количественной интерпретации.

Рис. 2.

Скоростной разрез земной коры в изолиниях скорости продольных волн V_P(x, z) по объединенному профилю Красноуральский + Ханты-Мансийский.

3. СЕТОЧНАЯ ПЛОТНОСТНАЯ МОДЕЛЬ ПО ХАНТЫ-МАНСИЙСКОМУ ПРОФИЛЮ

Наиболее подходящим для двумерного и трехмерного плотностного (гравитационного) моделирования является сеточный формат данных о распределении коровых скоростей по совокупности представленных разрезов ГСЗ. Границы слоев и блоков будущей плотностной модели выбираются по изолиниям сглаженной скоростной модели с тем, чтобы сохранить ее основные структурные особенности (рис. 2). Дальнейшие этапы вычислений демонстрируются на примере Ханты-Мансийского профиля длиной около 500 км. Обработка остальных четырех профилей (Красноуральский, Красноленинский, Тура–Орск и Сосьва–Ялуторовск) с построением двумерных плотностных моделей земной коры и верхней мантии выполнена в той же последовательности.

Скоростные разрезы даны в двумерном варианте. В связи с чем двумерная схема интерпретации профильных данных аномального гравитационного поля является вполне оправданной. Геометрия и размер сеток скоростных и плотностных матриц не изменяются. При заполнении плотностных ячеек уточняются только коэффициенты поинтервальной кусочно-линейной регрессии “скорость–плотность” $\sigma (V) = aV + b.$

Вычисления модельного гравитационного поля проводились с помощью двумерного интеграла Пуассона. Аномалии вертикальной составляющей гравитационного поля g_z(далее Δg) определяются с точностью до постоянной (или слабопеременной) величины g₀ – нормального поля равновесного-гидростатического сфероида с эквивалентным распределением масс. Соответственно, для 2D “плоских” плотностных моделей вместо абсолютной плотности $\sigma (x,z)$ будем рассматривать избыточную плотность $\Delta \sigma = \sigma (x,z) - {{\sigma }_{0}}(z)$, вычисляемую относительно гидростатической плотности ${{\sigma }_{0}}(z)$:

При этом следует обратить внимание на одну важную особенность задачи гравитационного моделирования. Плотностные модели земной коры и верхней мантии (верхней части литосферы) чаще всего имеют форму бесконечного по простиранию пласта с горизонтальными границами обрамления [Павленкова и др., 1991; Романюк, 1995; Куприенко и др., 2007; Дружинин и др., 2014]. Такое ограничение является непременным условием для вычисления гравитационных аномалий от неоднородного слоя с точностью до постоянной составляющей фона. Если плотностной разрез ограничен снизу криволинейной границей (границей М), то в расчетной аномалии поля появится переменная составляющая, амплитуда которой пропорциональна скачку плотности на криволинейной границе “кора–мантия”. Но вместо мантийной плотности в модель будет входить средневзвешенная плотность земной коры [Мартышко и др., 2010; Ладовский и др., 2017]. Поэтому нижнюю часть разреза (ниже границы М) следует дополнить горизонтальной границей, а во вновь образованном слое подобрать плотности мантийных масс.

При построении плотностного разреза по Ханты-Мансийскому профилю начальное распределение плотности $\sigma \left( {x,z} \right)$ выбиралось по двумерной скоростной матрице ${{V}_{P}}(x,y)$: выше границы М – на основе поинтервальной эмпирической зависимости “скорость–плотность”; ниже границы М плотность принимается постоянной с контролирующим ограничением ${{\sigma }_{M}} \geqslant {\text{3}}{\text{.3 г/с}}{{{\text{м}}}^{{\text{3}}}}$ при ${{V}_{P}} \geqslant \,\,{\text{8 }}\,{\text{км/с}}$ [Дружинин и др., 1982; 2014]. По плотностной матрице $\sigma \left( {x,z} \right)$ были вычислены значения средних для каждого гипсометрического уровня; найденные средние значения приравниваются к распределению зависящей только от глубины гидростатической плотности ${{\sigma }_{0}}\left( z \right)$. За пределами профиля плотностная модель дополнялась двумя бесконечными полупластами, плотность которых приравнена к гидростатической.

Далее, прямая задача гравиметрии (1) от двумерных плотностных моделей рассчитывалась по специально разработанной программе “V-σ-CALC” для скоростных сеточных матриц. По каждому скоростному интервалу кусочно-линейной зависимости “скорость–плотность” строилась матрица избыточных плотностей и вычислялся прямой гравитационный эффект модели до глубины 80 км – предполагаемого уровня региональной изостатической компенсации. При этом никакого априорного деления модели на слои и блоки не требуется, также не требуется вычисления фоновой составляющей поля от “законтурных масс”. Разрез, избыточные плотности которого вычисляются относительно гидростатической, находится в “пустом” пространстве, т.е. пространстве с нулевым гравитационным эффектом. Поле, рассчитанное по матрице избыточных плотностей, будет всегда приведено к нулевому уровню [Ладовский и др., 2017].

Расчетное поле консолидированной коры с однородной мантией, плотность которой приравнена к средневзвешенной плотности коры ${{\bar {\sigma }}_{K}}$, гораздо ближе к наблюденным аномалиям Буге. Если средневзвешенную плотность ниже границы М до глубины 80 км заменить мантийной (${{\sigma }_{M}} = {\text{3}}{\text{.3 г/с}}{{{\text{м}}}^{{\text{3}}}}$), то расхождение наблюденного и вычисленного полей становятся больше (рис. 3). Увеличение амплитуды длинноволновых аномалий модельного поля свидетельствует о возможности региональной компенсации мантийных плотностей ниже границы М. Граничное условие изостазии на глубине 80 км позволяет наметить контуры мантийных блоков изостатической компенсации и рассчитать для них значения плотности.

Рис. 3.

Расчетные гравитационные аномалии над Ханты-Мансийским разрезом с однородной мантией до глубины 80 км: 1 – проекция наблюденного поля на направление профиля; 2 – поле консолидированной земной коры без осадочного чехла и мантии ${\text{(}}{{\bar {\sigma }}_{K}} = {\text{2}}{\text{.87 г/с}}{{{\text{м}}}^{{\text{3}}}}{\text{)}}$; 3 – поле коры с осадочным чехлом и верхней мантией ${\text{(}}{{\sigma }_{M}} = {\text{3}}{\text{.3 г/с}}{{{\text{м}}}^{{\text{3}}}}{\text{)}}{\text{.}}$

4. МАНТИЙНЫЕ БЛОКИ

При построении мантийных блоков по матрице избыточной плотности создается сеточный файл распределения литостатического давления до заданной глубины. Литостатическое давление на глубине h равно весу неоднородной плотностной колонки $\sigma \left( {x,z} \right)$ единичного сечения; гидростатическое давление на той же глубине равно весу однородной колонки с гидростатической (одномерной) плотностью ${{\sigma }_{0}}\left( z \right)$. Отклонение $\Delta P\left( {x,h} \right)$ литостатического давления $P\left( {x,h} \right)$ от его среднего (гидростатического) значения ${{P}_{0}}\left( h \right)$ на глубине $h$ рассчитывается по формуле:

Аномалии литостатического давления пропорциональны избыточной плотности, так что плотностной разрез для модели с однородной мантией легко перестраивается в литостатический. В поле аномалий литостатического давления отчетливо конфигурируется блочная модель кристаллической земной коры и верхней мантии на разных глубинных срезах. Поправки, введенные в рамках предположения о частичной изостатической компенсации, позволяют перераспределить плотность в мантии таким образом, чтобы значение избыточного давления на предполагаемом уровне ${{H}_{{{\text{ГР}}}}} = 80\,\,~{\text{км}}$ стало близким к нулю: $~\Delta P\left( {x,{{H}_{{{\text{ГР}}}}}} \right) = 0.$ Для коррекции плотности в мантийных блоках была введена функция-компенсатор ${{\Sigma }}\left( x \right)$, которая показывает, какое значение плотности нужно добавить (или вычесть) в слое между границей М и уровнем ${{H}_{{{\text{ГР}}}}} = 80$ км, чтобы модель была бы изостатически скомпенсирована [Martyshko et al., 2017].

Пусть $\Delta {{P}_{{{\text{одн}}}}}$ и $\Delta {{\sigma }_{{{\text{одн}}}}}$ – отклонения литостатического давления и плотности от их усредненных (гидростатических) значений для модели с однородной мантией. Тогда при изменении плотности в мантии на ${{\Sigma }}\left( x \right)$ аномальное литостатическое давление $\Delta {{P}_{{{\text{бл}}}}}\left( {x,h} \right)$ в модели с блочной мантией на глубине изостатической компенсации h будет близко к нулю:

На рис. 4 приведена последовательная схема построения плотностной модели блочной мантии, согласно принципу изостазии. Скоростная модель с однородной мантией и ее плотностной эквивалент (рис. 4а) порождает столбообразные аномалии литостатического давления ниже границы М; на глубине предполагаемого уровня изостазии 80 км они меняются до 450 бар (рис. 4б). Функция-компенсатор (2) демонстрирует изменение по разрезу изостатических добавок к мантийной плотности (рис. 4в), распределение компенсирующей плотности ${{\Sigma }}\left( x \right)$ весьма неоднородно и по вариациям не соотносится с геофизическими моделями крупноблочного строения верхов мантии. Чтобы устранить это несоответствие, вертикально расслоенная модель ниже границы М была разбита на несколько блоков, внутри которых проведено усреднение плотности. X-координаты разбиения выбирались по нулям функции ${{\Sigma }}\left( x \right)$ с тем расчетом, чтобы получившиеся блоки не имели слишком малую протяженность вдоль профиля. При помощи усредненной таким образом ${{\Sigma }}\left( x \right)$ была построена блочная модель верхней мантии (рис. 4г). Для Ханты-Мансийского профиля латеральная протяженность мантийных блоков неплохо согласуется с интервалами сравнительно постоянных граничных скоростей головных волн на кровле верхней мантии [Дружинин, 1983].

Рис. 4.

Распределение аномалий литостатических нагрузок вдоль профиля Ханты-Мансийский: 1 – скоростная модель с однородной мантией (σ_М = 3.3 г/см³); 2 – блочное распределение аномалий литостатического давления $\Delta P\left( {x,h} \right)$; 3 – функция-компенсатор литостатических аномалий ниже границы М; 4 – модель с блочной мантией, построенная при условии изостатической компенсации масс на глубине 80 км. Цвет изолиний разреза земной коры соответствует скоростной и плотностной шкалам на рис. 2 и рис. 3.

5. ЛИНЕЙНАЯ ОБРАТНАЯ ЗАДАЧА ГРАВИМЕТРИИ

В процессе решения линейной обратной задачи гравиметрии подбираются не плотностные элементы сеточных матриц, число которых чуть меньше 4 × 10⁴, а только несколько коэффициентов кусочно-линейной регрессии и плотности мантийных блоков. Такое значимое уменьшение количества искомых неизвестных (параметров подбора плотностной модели) значительно повышает алгоритмическую устойчивость решения обратной задачи гравиметрии и сокращает затраты машинного времени. Это принципиальное отличие нашего подхода от тривиального подбора неизвестных плотностей на множестве разбиений элементарных тел сеточного каркаса.

Распределение скорректированных (подобранных) плотностей по скоростным разрезам находится из решения линейной двумерной обратной задачи гравиметрии. Одновременно подбираются плотности мантийных блоков:

Здесь: $g\left( x \right)$ – наблюденное поле; $F\left( {x,{{\sigma }_{1}},{{\sigma }_{2}}, \ldots ,{{\sigma }_{N}},{{\sigma }_{{{{B}_{1}}}}},{{\sigma }_{{{{B}_{2}}}}}, \ldots ,{{\sigma }_{{{{B}_{K}}}}}} \right)$ – модельное поле (1) в точке $x;~\,{{\sigma }_{{{{B}_{1}}}}}, \ldots ,{{\sigma }_{{{{B}_{K}}}}}$ – плотности $K~$ мантийных блоков. Переменные ${{\sigma }_{1}}, \ldots ,{{\sigma }_{N}}$ – это значения плотности, которые соответствуют скоростям продольных волн ${{V}_{1}}, \ldots ,{{V}_{N}}$. Они связаны с поинтервальными коэффициентами A_P и B_P посредством системы линейных уравнений:

Примечание: для Ханты-Мансийского разреза N = 4; $K = 5.$

Аномалии гравитационного поля пропорциональны плотностным разностям по заданным трем “P-скоростным” интервалам; плотности мантийных блоков отнесены к среднему значению плотности однородной мантии $\sigma _{M}^{0}~\,\, = \,~3.3\,\,~{{\text{г}} \mathord{\left/ {\vphantom {{\text{г}} {{\text{с}}{{{\text{м}}}^{{\text{3}}}}}}} \right. \kern-0em} {{\text{с}}{{{\text{м}}}^{{\text{3}}}}}}$:

Плотностные и скоростные разности связаны между собой только угловыми коэффициентами A_P линейной регрессии (4). Таким образом, только восемь неизвестных будут искомыми переменными в задаче минимизации (3): три коэффициента A_P и пять избыточных плотностей мантийных блоков. Этого достаточно, чтобы уточнить матрицу избыточных плотностей разреза расширенной сводной модели “кора–мантия”. Поинтервальные плотности разреза (в абсолютной мере) пересчитываются по разностям (5). Три коэффициента В_P находятся по разделенным значениям одномерной гидростатической плотности ${{\sigma }_{0}}(z)$ по “P-плотностным” интервалам:

Подобранные сеточные плотности по Ханты-Мансийскому разрезу до глубины 80 км сопоставляются с сеточным массивом исходных данных по скоростям продольных волн в коре и приграничной мантии. Связь двух однотипных числовых массивов позволила рассчитать поинтервальные значения коэффициентов материальных уравнений по результатам математического моделирования и уточнить исходную регрессионную зависимость “плотность–скорость” (скомпилированную по материалам ранее опубликованных работ [Дружинин и др., 1982] или полученную по петрофизическим выборкам на образцах горных пород [Мартышко и др., 2016а]):

На рис. 5 показана 2D-слоисто-блоковая сейсмоплотностная модель земной коры и верхней мантии по профилю Ханты-Мансийский. Модель построена по подобранной матрице избыточных плотностей, заверка коэффициентов формулы пересчета (6) выполнена по скоростной матрице исходных данных. Расчетное гравитационное поле модели сопоставлено с наблюденными аномалиями Буге. При этом постоянная составляющая фона g₀ не учитывается. Морфологическое сходство двух кривых для двумерных моделей вполне приемлемо на качественном уровне интерпретации.

Рис. 5.

Подобранный плотностной разрез по Ханты-Мансийскому сейсмическому профилю. Над разрезом приведены графики аномалий гравитационного поля, отнесенные к нулевому уровню: синяя кривая – фактические (наблюденные) данные, красная кривая – расчетные значения поля от модели неоднородного пласта до глубины 80 км. Условные обозначения тектонических структур: ВУПр – Восточно-Уральское погружение; ЗУП – Зауральское поднятие; ТКПр – Тюмено-Кустанайский прогиб; ХМСП – Ханты-Мансийское срединное поднятие; ЗСГС – Западно-Сибирская геосинеклиза (Фроловский блок).

Контуры мантийных блоков, намеченные по изостазии, отражают блочную делимость вышележащих структур верхней и средней земной коры. Интересно, что латеральная изменчивость плотности, конфигурирующая блочный разрез по Ханты-Мансийскому профилю, получена по результатам количественной интерпретации сейсмических данных и решения обратной задачи гравиметрии. Над разрезом (рис. 5) нанесены структуры из сводной тектонической схемы Урало-Сибирского региона (рис. 1). Отметим, что математическая и геологическая интерпретация практически согласуются, но количественная интерпретация позволяет также проследить структуры, выделенные по фундаменту, на глубину. Границы тектонических зон чередующихся поднятий и погружений консолидированного фундамента соответствуют блочной структуре плотностного разреза коры и мантии (рис. 5). Областям поднятий (Зауральское поднятие и Ханты-Мансийское срединное поднятие) соответствуют более легкие блоки консолидированной коры и более тяжелые – верхней мантии. В зонах погружений (Восточно-Уральский прогиб и Тюменско-Кустанайский прогиб) ситуация обратная: тяжелым блокам земной коры соответствуют более легкие блоки верхней мантии. Анализ форм рельефа границы фундамента и границы верхней мантии обнаруживает устойчивую тенденцию к образования антиформ кровли фундамента и кровли верхней мантии, т.е. обрамляющего внешнего контура консолидированной земной коры.

6. ТРЕХМЕРНЫЙ ФОРМАТ ПЛОТНОСТНОЙ МОДЕЛИ

Для построения трехмерной сейсмоплотностной модели Североуральского сегмента в пределах градусной трапеции (58°–64° с.ш., 54°–72° в.д.) были использованы данные сейсмических исследований, охватывающих большую площадь. Увеличение длины сейсмических профилей предполагает построение пространственного каркаса оцифрованных разрезов для расширенной территории с последующим замыканием контура интерполяции на требуемый сегмент. Это позволяет избежать появления высоко градиентных аномалий расчетного поля на границе планшета при сопутствующем разрыве продолженных плотностей во внешность трехмерной модели. По методике, изложенной в предыдущих разделах, в дополнение к Ханты-Мансийскому профилю были построены двумерные плотностные разрезы по четырем сейсмическим профилям: Вижай–Нижняя Тура–Орск, Северная Сосьва–Ялуторовск, Красноуральский и Красноленинский.

Использование двумерных плотностных моделей для перехода к трехмерной конфигурации не требует высокоточного подбора гравитационного поля по отдельно взятым профилям. Вполне достаточно качественного совпадения расчетного поля разреза с ортогональной проекцией на профиль трехмерных аномалий наблюденного поля. В процессе решения линейной обратной задачи гравиметрии подбираются не плотностные элементы сеточных матриц, размером $ \approx 1000 \times 800$ ячеек, а только угловые коэффициенты кусочно-линейных функций (4) по выделенным скоростным интервалам. Для Ханты-Мансийского профиля решение обратной задачи гравиметрии (3) построено для трех интервалов (6) изменения скорости в земной коре и пяти мантийных блоков. Для любого из оставшихся четырех плотностных разрезов количество искомых параметров обратной задачи меняется несущественно: число выбираемых скоростных интервалов в коре не более 3; число мантийных блоков (от 3 до 7) зависит от “нулей” функции компенсатора (2) изостатической редукции (см. рис. 4). Такое значимое уменьшение количества параметров подбора плотностной модели повышает алгоритмическую устойчивость решения обратной задачи гравиметрии и значительно снижает затраты машинного времени. В табл. 1 приведены уточненные коэффициенты кусочно-линейной зависимости “плотность–скорость” по каждому из профилей и мера несогласия (средне-квадратичное отклонение RMS в мГал) модельного и наблюденного полей. Оценки даны как для плотностных разрезов с однородной мантией RMS_од (начальное приближение), так и для подобранных разрезов с мантией блочной структуры RMS_бл.

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

Рис. 6.

Пространственное положение плотностных разрезов на координатном цифровом планшете карты аномалий гравитационного поля. Цифры в кружочках – профили Баженовской геофизической экспедиции: 1 – Красноленинский, 2 – Красноуральский и Ханты-Мансийский, 3 – Вижай–Нижняя Тура–Орск, 4 – Северная Сосьва–Ялуторовск.

Недостающие данные по плотностям в межпрофильном пространстве заполняются интерполированными значениями плотности с опорных разрезов трехмерного каркаса. Интерполяция выполняется по отдельным горизонтальным слоям. Мощности слоев определяются выбранным нами шагом дискретизации (100 м) по оси глубин; в горизонтальных плоскостях шаг интерполяции 1 км. Результатом интерполяции является четырехмерный массив послойных сеточных файлов (x, y, z, σ), в котором каждому горизонтальному слою сопоставлена глубина, на которой он расположен. Данные за пределами контура интерполяции не учитываются. Результатом работы программы интерполяции является пакет сеточных “грид-файлов” (послойных горизонтальных сечений $z = {\text{const}}$) четырехмерного (x, y, z, σ) массива и XML-файл, в котором каждому слою сопоставлена глубина, на которой он расположен:

Построенный “цифровой параллелепипед” на рис. 7 представляет собой трехмерную плотностную модель начального приближения для заданного участка площади. Алгоритм “точного интерполятора” настраивается на трехмерную конфигурацию (рис. 6): он учитывает схему тектонического районирования и, по возможности, зональность гравитационного поля. Непрерывное послойное заполнение сеточных матриц “цифрового параллелепипеда” восстанавливает интерполированные значения плотности в межпрофильном пространстве и сохраняет подобранные 2D-плотностные параметры в сечении сейсмических профилей. Градиентно-слоистая модель интерполированной плотности может рассматриваться как первое (начальное) приближение трехмерной модели земной коры и верхней мантии. Модель легко перестраивается для учета дополнительной информации и доступна для редактирования с возможностью получения информации по любому сечению вдоль образующих кривых, заданных в специальном формате граничного файла “.bln”. На основе послойных сеточных матриц трехмерного параллелепипеда реализуется алгоритм построения карт горизонтальных срезов и структурных карт рельефа опорных плотностных границ, что принципиально важно при решении структурных задач тектонического районирования.

Рис. 7.

Плотностная 3D-модель начального приближения, полученная интерполяцией плотностей с двумерных разрезов в межпрофильное пространство. Сверху нанесена карта расположения сейсмических профилей и схема тектонического районирования (см. рис. 1).

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

7. ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

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

Практические задачи, возникающие при 3D-моделировании глубинного строения неоднородной геологической среды, требуют минимизации ручного труда не только на этапе подготовки исходных данных, но и на всех этапах промежуточных расчетов. Формализация исходных данных ГСЗ и их преобразование в сеточный формат скоростных (и плотностных) разрезов способствует разработке высокоэффективных вычислительных алгоритмов “быстрого” решения прямых и обратных задач гравитационного моделирования на сетках большой размерности [Мартышко и др., 2013; Martyshko et al., 2018]. Доминантой предполагаемых решений служит плотностная модель начального приближения. По этой причине мы скрупулезно выполнили все промежуточные вычисления, начиная от томографии ГСЗ (поля времен и сеточные модели скоростных разрезов) и заканчивая трехмерной моделью интерполированной плотности. Кратко перечислим необходимые промежуточные этапы.

1. По специальным двумерным полям времен рефрагированных продольных волн в первых вступлениях строится томографический скоростной разрез земной коры.

2. По сеточной скоростной матрице и кусочно-линейной петрофизической зависимости “скорость–плотность” рассчитывается массив предварительных значений плотности на сетке той же размерности.

3. Модель земной коры дополняется однородным мантийным слоем с плотностью ${{\sigma }_{M}} = {\text{3}}{\text{.3 г/с}}{{{\text{м}}}^{{\text{3}}}}$ до глубины регионального уровня изостатической компенсации 80 км.

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

5. Строится горизонтально-слоистая (одномерная) модель “нормального” плотностного разреза. Одномерная, усредненная по горизонтальным сечениям, плотность ${{\sigma }_{0}}\left( z \right)$ принимается за гидростатическую.

6. Плотностной разрез трансформируется в модель горизонтального пласта бесконечного простирания. Для этого внешнее законтурное пространство за пределами профиля заполняется массами с одномерной гидростатической плотностью.

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

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

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

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

8. ЗАКЛЮЧЕНИЕ

Предложена технология преобразования скоростных разрезов в двумерные плотностные матрицы. Разработан алгоритм трансформации двумерных числовых массивов в формат послойных сеточных функций трехмерного плотностного куба. Построены сейсмоплотностные разрезы вдоль пяти профилей ГСЗ Североуральского сегмента. На сейсмических разрезах выделены опорные границы, характеризующие поверхность фундамента и подошву земной коры (границу Мохоровичича), а также внутрикоровые поверхности. Уже на первоначальном этапе анализа построенных двумерных моделей просматриваются общие и отличительные признаки распределения взаимоувязанных плотностных и скоростных параметров по глубине. В верхней и средних частях земной коры уверенно трассируются зоны сочленения Северного Урала с Восточно-Европейской платформой и Западно-Сибирской плитой. Значимые оценки получены и для усредненной плотности мантийных блоков в пределах указанных структур. Таким образом, разработанная методика двумерного сейсмоплотностного моделирования позволяет учесть основные структурные особенности модели скоростного разреза ГСЗ и реализовать устойчивый подбор плотностей в классе элементарных тел (элементов сеточной матрицы) для любой, сколь угодно сложной, геометрии среды. Конвертация разнородных числовых массивов к стандартному (матричному) формату двумерных и трехмерных сеточных функций способствует созданию новых эффективных алгоритмов математического сопровождения и рационализации вычислительных процессов. Все разработанные алгоритмы реализованы в прикладных программах с использованием параллельных алгоритмов для многопоточных вычислений. По пяти профилям (Ханты-Мансийский, Красноуральский, Вижай–Нижняя Тура–Орск, Северная Сосьва–Ялуторовск и Красноленинский) построены двумерные модели распределения плотности, для которых выбором конфигурации мантийных блоков удалось с хорошей точностью выполнить условие изостатического равновесия на глубине 80 км. Полученные разрезы в дальнейшем будут использованы для построения модели объемного распределения плотности на трехмерных сетках большой размерности и позволят получить новые структурные карты поверхностей фундамента и Мохо для Североуральского сегмента.