Обратная задача моделирования тектонической эволюции рифтовых бассейнов тема автореферата и диссертации по физике, 01.04.12 ВАК РФ

?Г6 ОА - 8 ОМ

^ЩЩОНАЛЬНАЯ АКАДЕМИЯ НАУК УКРАИНЫ Институт геофизики им. С.И. Субботина

На правах рукописи УДК (550.36+551.24.03):550.83.017

Поплавский Константин Николаевич

ОБРАТНАЯ ЗАДАЧА МОДЕЛИРОВАНИЯ ТЕКТОНИЧЕСКОЙ ЭВОЛЮЦИИ РИФТОВЫХ БАССЕЙНОВ

лес

I I Г —- (

Специальность: 01.04Л2 - Геофизика

Автореферат диссертации на соискание ученой степени кандидата физико-математических наук

Киев - 1996

Диссертация является рукописью. Работа выполнена в Институте геофизики им.С.И.Субботина НАЛ Украины

Научный руководитель:

Официальные оппоненты:

академик HAH Украины,

доктор физико-математических наук,

профессор В. И. Старостенко

доктор физико-математических наук, профессор Г.Я. Голиздра

кандидат физико-математических наук Я.М. Хазан

Ведущая организация:

Днепропетровское отделение Украинского государственного Института минеральных ресурсов, г. Днепропетровск

Защита состоится "_" _ 1996 г. в_часов

на заседании специализированного совета Д 01.95.01 при Институте геофизики им.С.И.Субботина НАН Украины: 252680, г. Киев-142, пр. Палладина, 32.

С диссертацией можно ознакомиться в библиотеке Института геофизики им.С.И.Субботина HAH Украины.

Автореферат разослан "_" _ 1996 г.

Отзывы на автореферат просим направить в двух экземплярах ученому секретарю специализированного совета по адресу: 252680, г. Киев-142, пр. Палладина, 32.

Ученый секретарь

специализированного совета У .^-"у

доктор физико-математических наук II ллУ" / В.С Гейко.

[ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Ашщдышат_работы: Рифтовые бассейны составляют обширный и важный подкласс осадочных бассейнов. Многие из них содержат значительные запасы углеводородов. Численное моделирование механизмов возникновения рифтовых бассейнов и их развития представляет значительный исследовательский и практический интерес.

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

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

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

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

* сконструирована общая схема решения обратной задачи моделирования;

* реализована и опробована итерационная схема одновременной многопараметрической минимизации глобальной целевой функции;

* выделен и формализован (с использованием принципа коллокации по областям) набор взвешенных целевых функций на основе анализа физического смысла поставленной задачи;

* разработаны и реализованы эффективные алгоритмы индивидуальной минимизации взвешенных целевых функций;

* разработан и реализован технологичный программный комплекс, пригодный для решения как исследовательских, так и практических задач;

* разработанный подход опробован при моделировании реального геологического объекта.

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

Прдтшнеское^начшие м^тжм1ищ^ЩАШ1шат-исследшш1шй. При исследовательском моделировании предлагаемый способ предоставляет возможность проанализировать влияние вариаций многочисленных входных параметров алгоритма прямой задачи. При этом исследователь легко сохраняет контроль над реализациями модели. Это превращает предлагаемый способ в эффективный инструмент исследования как реализаций модели, так и самих модельных представлений, позволяя анализировать пределы эквивалентности решений. Эффективность и управляемость открывают возможность дальнейшего развития алгоритмов прямой задачи в направлении подключения тектонических процессов, ранее игнорировавшихся. При ручном моделировании этот путь развития был практически исчерпан, так как оператор прямой задачи с ростом числа параметров быстро становится неконтролируемым для интерпретатора.

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

Оснйвныетщищае.мыелаучныеположенит

1) Разработанный и реализованный способ решения обратной задачи моделирования тектонической эволюции рифтовых бассейнов позволяет более эффективно и надежно выполнять реализацию модельных представлений интерпретатора.

2) Защищаемый способ предоставляет интерпретатору эффективный инструмент исследования модельных представлений.

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

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

Апробация работы: Основные положения работы докидывались на следующих конференциях и семинарах:

1. EUROPROBE Workshop, Ворзель, Украина, октябрь 1994 г.

2. EAUG-8 annual meeting, Strasbourg, France, April 1995.

3. EUROPROBE Workshop, Leeds, United Kingdom, My 1995.

4. AUG 1995 fall meeting, San Francisco, USA, December 1995.

5. Structural geology and tectonics seminar, Institute of Earth Sciences, Vrije Universiteit, Amsterdam, The Netherlands, September 1995.

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

- в исследованиях, проводимых в отделе глубинных процессов Земли и гравиметрии Института геофизики им. С. И. Субботина HAH Украины.

- в работах, проводимых на кафедре тектоники и структурной геологии Свободного Университета (Vrije Universiteit, Amsterdam, The Netherlands).

Публикации! По теме диссертации опубликовано 7 работ.

Объем и структура работы: Диссертация состоит из ВВЕДЕНИЯ, глав СОСТОЯНИЕ ПРОБЛЕМЫ И ФОРМУЛИРОВКА ЗАДАЧИ, СПОСОБ РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ МОДЕЛИРОВАНИЯ, ПРИМЕР И РЕЗУЛЬТАТЫ ПРИМЕНЕНИЯ, ЗАКЛЮЧЕНИЯ и ПРИЛОЖЕНИЯ, общим объемом 132 страницы. Работа включает 2 таблицы, 32 рисунка, содержит список литературы из 112 наименований.

Фактический материал и личный вклад автора: В основу диссертации легли исследования, выполнявшиеся автором в Институте геофизики им. С.И.Субботина HAH Украины а также в Свободном Университете (Vrije Universiteit, Amsterdam, The Netherlands) на кафедре тектоники и структурной геологии. Использовавшиеся при моделировании наблюденные данные (сейсмостратиграфические разрезы по профилям в рамках

-3-

международного проекта EUROPROBE) являются результатом многолетнего труда коллектива специалистов ГГП "Укргеофизика".

Работа выполнялась под руководством академика HAH Украины, доктора физико-математических наук, профессора В. И. Старостенко, которому автор выражает глубокую и искреннюю благодарность.

Автор выражает большую признательность к.ф.-м.н. Ю.Подладчикову, без деятельного участия которого эта работа не была бы ни начата, ни закончена. Неоценимую пользу автору принесли обсуждения, консультации и советы доктора Р.Стифенсона (Dr. R.A. Stephenson), доктора Х.Куи (Dr. H.Kooi), д-ра г.-м. наук В.Г. Гутермана, чл.-кор. НАНУ ВАДаниленко, Д.Б.Венгровича, профессора Н.Кузнира (Prof. N.J.Kusznir), профессора С.Клутинга (Prof. S.A.P.L.Cloetingh), чл.-кор. НАНУ П.Ф.Шпака, профессора Д.Юна (Prof. D.A.Yuen), к.г.-м.н. С.Н.Стовбы и других коллег, которым автор сердечно благодарен.

КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ

^ВВЕДЕНИЕ

Во ВВЕДЕНИИ рассмотрены актуальность изучаемой проблемы, цели," задачи, научная новизна и практическое значение работы. Изложены основные защищаемые положения.

2.1 Краткий обзор развития численных моделей эволюции рифта Приводится разделение концепций рифтинга на гипотезы "активного" и "пассивного" рифгогенеза. Основная задача обзора - представить историю формирования блоковой и изгибно-консольной моделей. Ограничение этими двумя моделями связано с тем, что только они (из известных автору) ориентированы на двумерное моделирование осадочного заполнения бассейна, т.е. стратиграфии осадков, а именно способ анализа по-стрифтовой стратиграфии и является целью данной работы.

2.2 Аиоршпмы решения прямой задачи моделирования рифтинга

Раздел содержит описание трех алгоритмов решения прямой задачи моделирования эволюции рифта - модели МакКензи, блоковой модели и изгибно-консольной модели. Модель МакКензи является исторически первой и базовой для всего описываемого подхода. Краткое описание изгибно-консольной модели приводится для сравнения.

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

- 4 -

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

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

а

/ 3 2,.Л 32

дх

d2w

л(х'0ч~т + N Д-Г+ (р„ - pf)gw(x) = q(x,t);

, л , -1-1 IJ — I/ , I1-1VI Д I - UI I. it. ^ ^

dx, I дх'

Здесь и далее обозначения, не объясненные в тексте, см. в Таб. 1. Изгибная жесткость литосферы D определяется как изгибная жесткость фиктивной упругой плиты мощностью ЕЕТ. , . _ Е •(EET(x,t))

{X,t>~ 12(1-v2) ЕЕТ в каждом блоке определяется глубиной до изотермы, задаваемой интерпретатором. Суммарная вертикальная нагрузка определяется как:

<?(*)= J (й a (z)~ p(*,Z,t ))dz I i2'2)

о

где р0(z) - начальное распеределение плотности в нерастянутой литосфере, а p(x,z,t)- распределение плотности в литосфере и слое осадочного заполнения депрессии после растяжения, остывания в течение времени tk и заполнения начального погружения осадками.

Изменение плотности вещества в результате контракции и положение изотермы, определяющей ЕЕТ, вычисляется после нахождения распределения температур через заданное время tk при решении двумерного уравнения теплопроводности (методом конечных разностей):

д2Г д2Т +

дх2 dz

2

( 2-3 )

у

дТ

Граничными условиями при решении этого уравнения устанавливаются значения То-0 при 1~0 и Т=Та при г=/г£ . Начальные условия при старте алгоритма определяются геометрией модели и заданными в каждом блоке величинами растяжения. Теплогенерация в литосфере, объемная контракция и теплоизолирующее влияние осадков игнорируются. Го- 5 -

ризонтальные нагрузки задаются интерпретатором в виде величин Л^ горизонтальных палеонапряжений.

Этапами работы алгоритма являются к возрастов стратиграфических слоев, имеющихся в наблюденном разрезе. Задача расчета на каждом этапе - определить мощность осадочного слоя, накапливающегося за к-с время, равную полному погружению литосферы за это время. Полное погружение 5(х) рассчитывается как сумма тектонического погружения 5ис(х) и упругой реакции литосферы на "догрузку" осадками иу (х). Тектоническое погружение Я,ес(х) рассчитывается как результат упругой изостатической компенсации тектонических нагрузок при нулевой плотности слоя осадочного заполнения. Для расчета упругого отклика литосферы на изменяющиеся величины осадочных нагрузок внутри каждого к-го этапа организованы вложенные итерации.

Рассчитанное таким образом погружение на к-ом этапе дает значение стратиграфической мощности к-го горизонта в /-том блоке при его отложении - обозначим эту мощность как . На следующем, к+1-м этапе в значения мощностей вводятся поправки за литофикацию осадков.

После последнего этапа расчета весь расчетный стратиграфический разрез известен - т.е. известны мощности каждого к-го горизонта в каждой у'-ой точке профиля для момента времени РВ - сегодняшнего дня. Пусть все эти мощности образуют вектор Г (вектор расчетной стратиграфии) длиной Му. Пусть также все входные параметры вышеописанного алгоритма прямой задачи (приведенные в Таб. 2 Приложения) образуют вектор X (вектор параметров модели) длиной Ых. Обозначим вышеописанную процедуру расчета стратиграфического разреза как ¥ и будем называть ее оператором прямой задачи. Оператор прямой задачи выполняется над вектором параметров модели. Его результатом является вектор расчетной стратиграфии (иначе, вектор состояния модели), то есть

2.3 Формулировка геофизической задачи

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

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

-6-

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

.?. / Общее описание подхода

Способ решения обратной задачи включает в себя алгоритм решения прямой задачи. Пусть С(Х) = У0-У=У0-Р(Х); тогда решению подлежит уравнение:

Где Р - оператор прямой задачи; X - вектор параметров, используемых применяемым алгоритмом решения прямой задачи; У - вектор состояния модели (вектор результатов, описывающий расчетный стратиграфический разрез бассейна); У0 - вектор наблюденного состояния модели (вектор наблюдений). Под наблюдениями при этом понимается наблюденный (на сегодняшний день) стратиграфический разрез вдоль моделируемого профиля. Вектор У0 состоит из наблюденных мощностей соответствующих стратиграфических горизонтов и поэтому в дальнейшем часто сокращенно упоминается как "наблюденная стратиграфия". Вектор Г, соответственно, в дальнейшем сокращенно упоминается как "рассчитанная стратиграфия".

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

Предположим, что точки У, являющиеся отображениями (посредством оператора прямой задачи) точек X области Ох пространства параметров, образуют в пространстве состояний модели область Оу Примем упрощающие предположения об однозначности и непрерывности отображения из области Ох в область Оу\ и о дифференцируемости функции У=Р(Х) по всем ее аргументам в области Оу.

Оба вектора ¥0 и К представляют собой точки в пространстве состояний модели. Мерой близости точек считается расстояние между ними, определяемое в метрике следующим соотношением:

Назовем такое расстояние "глобальной целевой функцией" или "глобальной ошибкой". Решением (3-4) считается любой вектор X, для которого расстояние отточки У=Р(Х) до точки Уа меньше заданного.

3.2 Постановка задачи

Дано:. Непрерывная область Ох допустимых значений параметров (элементов вектора X) в пространстве параметров; оператор F отображе-

(3-4)

(3-5)

ния вектора Хиз области Ох в область Оу пространства состояний модели (оператор прямой задачи); вектор У0 наблюденного состояния модели. Область Оу полагается областью дифференцируемости функции >'

Нужно найти: вектор Хг, принадлежащий области Ох , для которого расстояние между Г0 и Уг—Р(Хг) меньше заданной величины е.

Решаемое уравнение (3-4) является векторной записью системы из Иу уравнений с к* переменными. Будем искать его решение методом Ньютона с помощью итерационного процесса вида:

= (з-б)

Полагая, что + получаем, что приближение к решению на

каждой итерации заключается в построении переопределенной (при Л'у >Л'Х) системы линейных уравнений: С'(Х,)С,=-С{Х)- (3-7)

и ее решению относительно С1. Формирование матрицы Якоби (7 '(Х-,) производится с использованием аппроксимации частных производных

конечными разностями вида (поскольку (X ) = (X ) Цн*,) - {хи..., ХяЛ,Хл + А*., ^... ,ХИ,) - /и (*,)}

где т = 1, Nу; и ДХ - вектор приращений, задаваемый априорно. Итерационный процесс (3-6) реализуется итерационным алгоритмом с апостериорным заданием числа итераций. Необходимое для запуска итерационного алгоритма начальное приближение модели (начальный вектор Х0) строится на основе консервативного принципа наименьшей сложности модели. Наипростейшей моделью в данном случае является модель МакКензи. Начальное приближение для факторов растяжения строится из соотношений баланса отдельных колонок (блоков) в предположении локальной изостатической компенсации:

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

3,3 Опыт многопараметрической оптимизации

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

Для исследования алгоритма решения обратной задачи был сгенерирован синтетический разрез К , в вектор параметров X, которого затем были искусственно внесены искажения. Эти искажения должны были быть обнаружены и исправлены алгоритмом решения обратной

задачи. Искажения вносились в профиль фактора 8: профиль фактора Р; глубину некинга литосферы значения внутриплитных палеонапря-

жений Мк.

Процедуры одновременной многопараметрической оптимизации были построены на основе итерационного процесса (З-б).

Использовался и механизм штрафных функций в форме притягивания решения к "правильному" ответу. Тем не менее даже в таком предельно упрощенном случае не удалось добиться устойчивой сходимости алгоритма. Алгоритм быстро (в пределах нескольких итераций) сходится к окрестности минимума, а затем взрывообразно расходится. Были испробованы несколько других приемов, в частности, некоторые модификации метода скорейшего спуска, уменьшение размерности вектора параметров X, различные реализации штрафных функций и т.д.. Основные характеристики поведения алгоритма не менялись. Возможно, полезным оказалось бы усиление адаптивных свойств алгоритма в виде включения адаптивных инкрементов. Однако вычислительные затраты (как и при использовании методов статистической оптимизации) сделали бы алгоритм непригодным для практического моделирования.

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

3.4 Выбор целевых функций и способов оптимизации

Для каждого параметра (элемента вектора X ) будет выбираться своя, индивидуальная целевая функция и свой способ ее минимизации. Весовые коэффициенты при элементах глобальной целевой функции будут принимать (в соответствии с принципом коллокации) значения 1 или 0, в зависимости от того, входит данный элемент глобальной целевой функции в рассматриваемую частную (взвешенную) целевую функцию, или же нет.

Целесообразно при этом уменьшить размерность оптимизируемой части вектора X. Теперь в алгоритме решения обратной задачи вектор параметров X содержит не все параметры, используемые алгоритмом прямой задачи, а только параметры, представляющие первоочередной интерес для исследования реализаций модели (и наиболее трудоемкие для ручного подбора). К таковым были отнесены: а) профиль фактора 5; б) профиль фактора Р; в) глубина некинга литосферы 2песК; г) значения горизонтальных (внутриплитных) палеонапряжений Лгь действовавших во время накопления пострифтовых стратиграфических слоев.

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

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

Фактор 5 обеспечивает главный вклад в начальное погружение, определяющее, в свою очередь, конечную глубину до фундамента (полное погружение). Фактор (5 (точнее, его уклонение от фактора 8) определяет положение да кривой точки, соответствующей концу синрифтового периода. Наличие ненулевых пострифтовых палеонапряжений (1) приводит к ундуляции пострифтовой части кривой вокруг ее "нормального" экспоненциального убывания. Влияние изменений глубины некинга 2„еск будет рассмотрено ниже.

Фактору 6(х) ставится в соответствие глубина до фундамента (полное погружение, сумма всех мощностей): фактору Р(х) - глубина до подошвы пострифтовых осадков (пострифтовое погружение, сумма мощностей пострифтовых горизонтов); величине палеонапряжения Nк(1) для времени X - индивидуальная мощность отдельного пострифтового горизонта номер к, откладывавшегося во время !.

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

Целевые функции и оптимизация факторов 6 и р. В соответствии с анализом кривой погружения в качестве целевой функции для оптимизации параметра ё(х) выбирается расхождение между расчетными и наблюденными величинами полного погружения:

Практически это означает, что ищутся такие значения 8(х) , которые дают наилучшее совпадение расчетных и наблюденных глубин до фундамента, и только, не заботясь об остатьных элементах вектора ¥[.

(3-9)

В качестве целевой функции для оптимизации параметра ¡3(х) принимается расхождение между расчетными и наблюденными величинами пострифтового погружения:

Для одновременной оптимизации параметров 5(х) и ?>(х) был построен итерационный алгоритм, основанный на (3-6). Для вычисления частной производной по 5 фиксировались значения , а все значения 3} инкрементировались на одну и ту же величину Д8 . Значения самих частных производных д' по вычислялись отдельно для каждого блока:

^—^{г^у-зг'^^^лб^}: ( 3-11)

где ( Б ( 8,) ) - оператор извлечения полного погружения для /-того блока на /-той итерации. Аналогичная процедура применялась для вычисления частной производной по фактору р .

Полученные значения производных и величины приращений Д5 и др затем использовались для обновления

значений 5 и р на следующем шаге

итераций:

5,. ,, = 5. у + С,- г А8 = 5,) + 45 ( 3"12 )

Описанный итерационный алгоритм демонстрирует быструю и устойчивую сходимость не только по целевым функциям 5 и р факторов, но также по значениям глобальной ошибки. На надежную работу этого алгоритма во многом опираются последующие построения. Определим этот алгоритм как "оптимизирующий оператор" I (X, ¥{). При этом X'; —I1 ( У,) =11 (Г(X,)); где вектор X/ содержит оптимизированные профили 8(х) и р(х), в то время как все остальные параметры в нем остаются неизмененными. Оптимизирующий оператор I' позволяет построить итерационные схемы с вложенными итерациями, ре-оптимизирующими параметры первого порядка после внесения изменений в значения других параметров. Ре-оптимизация необходима для выявления ошибок решения, обусловленных вариациями параметров второго порядка, которые иначе будут экранированы гораздо более значимыми последствиями неоптимальности параметров первого порядка.

Итерационный алгоритм на основе оператора I1 устойчиво сходится к некоторому решению, которое, однако не является единственным по крайней мере по значениям р(х). Алгоритм прямой задачи, использующий двумерное уравнение теплопроводности, оказывается малочувствительным к горизонтальным вариациям фактора р и весьма чувствительным к его интегральному значению. Возникает потребность в регуляризирующем (в смысле ограничения класса допустимых решений) подходе.

В качестве такой процедуры используется фильтрующий оператор Ф ("оператор теплового фильтра"). Результат работы теплового фильтра

представляет собой профиль Р(х), приближенный по форме к функции сой(2х) с неизмененным интегральным значением Д Функция соь(2х) выбрана потому, что она является собственной функцией уравнения теплопроводности и дает распределение температур, рассасывающееся с наименьшей скоростью. Отфачьтрованный профиль эквивалентен по его атаянию на последующую эволюцию модели исходному (предложенному оптимизирующим оператором I1). В работе [11] приведено строгое обоснование эквивалентности. Интересно отметить, что глобальная ошибка нечувствительна к применению оператора Ф, а целевая функция фактора Р даже показывает убывание, т.е. применение теплового фильтра улучшает решение по крайней мере по ¡5.

Цшвая футщш^гутиызенита. Глубина некинга литосферы является важным параметром, изменение которого влияет на все мощности в расчетной стратиграфии, т.е. на все элементы вектора состояния модели К Выбор целевой функции для этого параметра не очевиден. Был проанализирован набор из четырех возможных функций-кандидатов (глобальной ошибки, целевой функции 5, целевой функции ¡3 и интегрального уклонения ¡5 от 5) для трех моделей - синтетической, реальной модели ЬК и реальной модели БЬ. Интегральное уклонение фактора /3 от 5 определяется как:

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

Из результатов анализа вытекает, что применение теплового фильтра по крайней мере не ухудшает качество получаемых решений. Очевидно также, что целевая функция 8 практически не чувствительна к глубине некинга в довольно широком диапазоне значений (примерно от 10 до 30 км). Следует подчеркнуть, что "нечувствительность" здесь означает, что при любом значении ^песк из этого интервала оптимизирующий оператор 1! способен добиться примерно одинакового минимума целевой функции 6, но профили 8(х) и/или р(х) при этом каждый раз будут разные. Из анализа поведения функций-кандидатов с привлечением критерия минимального усложнения модели вытекает, что линейная комбинация целевой функции ¡3, глобальной ошибки и уклонения Д от 6 представляет собой хорошую кандидатуру на роль целевой функции для параметра 2яеск. При этом выбор меньшей глубины некинга позволяет удержать модель ближе к простейшей возможной (конечно, внутри диапазона эквивалентности, т.е. не поступаясь при этом достигнутым минимумом расхождения между расчетной и наблюденной стратиграфиями).

( 3-13 )

Минимум линейной комбинации указал искомые оптимальные значения а именно: 10 км для синтетической модели, 10 км для модели ЬК и 17 км для модели 5Ь.

Целевая функция вмтин па-хтшпряжетй* В качестве целевой функции для оптимизации величины палеонапряжений N^,(0 , действовавших во время образования горизонта номер к , принимается расхождение между расчетными и наблюденными мощностями этого горизонта:

ль М

где - оператор изатечения мощности к -го горизонта в у -том блоке.

Величины палеонапряжений для каждого к-го возраста определяются регулярным перебором значений внутри предела прочности (задаваемого интерпретатором). Значение, минимизирующее целевую функцию палеонапряжений, затем принимается в качестве решения. Основным упрощающим предположением, позволяющим применить метод перебора, является предположение о том, что оптимизированное значение (предшествующий горизонт) не зависит от результата оптимизации значения Л^ (последующий горизонт). В противном случае перебор даже с грубым шагом потребовал бы совершенно неприемлемых вычислительных затрат. К счастью, это упрощающее предположение в комбинации с вышеуказанной целевой функцией работают вполне удовлетворительно. Более того, в этом случае, вероятно, может быть доказана теорема о сходимости (благодаря поиску решения на компакте набора параметров Ох и гладкому поведению функции упругого отклика литосферы на нагрузки).

3.5 Алгоритм решения обратной задачи моделирования Алгоритм решения обратной задачи моделирования организован с разделением на блоки: а) генерация начального приближения; б) однократный запуск оптимизирующего оператора /'; в) оценка диапазона эквивалентности глубины некинга; г) оптимизация величин палеонапряжений.

41_ЛРИМЕ£Л1РЕЗУЛЬТ/\1Ы^Г1РИ1У1Е11ЕШ1Я

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

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

ществующих взглядов на строение и историю формирования ДЦВ выходит за рамки работы. Поэтому содержание этого раздела ограничено сведениями обзорного характера, в минимальном объеме сопровождаемыми необходимыми обсуждениями.

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

4.2 Использованные профили и задана моделирования ДДВ занимает едва ли не уникальное положение среди мировых рифтовых структур в силу ряда обстоятельств. Среди них - высокая степень изученности (в т.ч. региональными профилями МОГТ). Два из этих профилей - Лосиновка-Кинашевка (1ЛС) и Сагайдак-Лебедин (БЬ) и были использованы для моделирования в настоящей работе.

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

Первые опыты двумерного численного моделирования стратиграфии ДДВ (в работе Кшгшг е1 а1,, 1996) показали, что с точки зрения "нормальной" эволюции рифта в ДДВ были аномально высокие скорости осадконакопления в карбоне, в начальную фазу пострифтовой истории. Для объяснения этого в упомянутой работе была привлечена концепция "мантийного плюма".

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

- 14 -

ние и развитое мантийного пяюма ("горячей точки", "астеносферного" или "мантийного диапира" и т.п.), который по определению является возмущением трехмерного теплового поля, с таким соотношением размеров. Равноправность пространственных осей в процессе развития диапира в XY-изотропной среде неизбежно должна привести к изометричности его формы в плане. Вопрос сводится к двум моментам: а) либо линейность рифта является свойством, привносимым в процесс взаимодействия литосферой, а не тектонической первопричиной; б) либо линейность рифта является результатом локализации деформаций, возникающих из-за напряжений регионального масштаба, изначально обладающих XY-анизотропией.

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

Если образование рифта как региональной структуры происходит под влиянием напряженно-деформированного состояния литосферы, то напряженное состояние литосферы должно оказывать воздействие я на его последующую (пострифтовую) эволюцию. Вполне логичным поэтому представляется анализ возможности объяснения искажений "нормальности" пострифговой кривой погружения за счет действия поля горизонтальных напряжений.

Целью анализа является исследование возможности объяснения увеличенных мощностей карбона в ДЦВ действием горизонтальных пале-онапряжений во время накопления пострифтовых осадков. Средством анализа является предлагаемый способ решения обратной задачи двумерного моделирования рифинга. Результаты предсташхены в следующих разделах.

4.3 Обсуждение результатов моделирования

Моделирование было выполнено для двух профилей (LK и SL) и заключалось в применении разработанной программы решения обратной задачи для нахождения таких значений параметров (профиль фактора S, профиль фактора р, глубина некинга литосферы, величины палеонапряжений), которые минимизировали бы расхождения между наблюденной и расчетной стратиграфиями вдоль моделируемого профиля. Другие параметры (см. Таб. 2 Приложения) были заданы главным образом из литературы. Результатом моделирования являются предложенные алгоритмом решения обратной задачи (оптимизированные) значения параметров, которые и подлежат последующей интерпретации.

Выбор^глуошшжкинга^ттасфлры, Консервативный принцип неусложнения модели диктует выбор наименьшего значения ' из диапазона эквивалентности (см. раздел 3.4). Кроме верности Оккаму, такой выбор имеет и содержательное преимущество: возможность задавать меньшие значения фактора ß, требуемые для получения стратиграфии, близкой к наблюденной. Полученные значения фактора ß вполне достаточны для объяснения интенсивных вулканических процессов и привлечения дополнительных причин не требуется. Имеется еще одно следствие: при выборе малой глубины некинга тем самым отдается предпочтение модели эволюции ДЦВ с малым либо вовсе отсутствовавшим син-рифтовым подъемом флангов. Существенно, что это предсказание модели является потенциально проверяемым (методом FTA) и, возможно, будет проверено в ближайшее время.

А/шаз _ра<яетньЕСЖж5лшкшшлеометрий_ бассейна, Сравнение наблюденной и расчетной стратиграфии показывает, что имеющиеся расхождения носят главным образом локальный характер. Эти расхождения не могут быть удовлетворены изменениями в параметрах компакции осадков (несмотря на визуальное впечатление), что проверялось численными экспериментами с использованием специальной реализации алгоритма решения обратной задачи. Наиболее вероятной причиной расхождений следует считать нарушение условия двумерности галокинезом. Расхождения в целом невелики, что позволяет считать задачу численного моделирования решенной.

Профили факторов растяжения имеют средние значения для профиля LK 5int=1.231 и ßint=2.997; для профиля SL 5int=1.393 и ß1Tlt =3.995. Отношение факторов 5int для профилей SL и LK, равное 1.13, весьма близко к отношению ширин центрального грабена ДЦВ в сечениях этих профилей, равному 1.14. Это представляется хорошим аргументом в пользу утверждения о преимущественно разломном режиме деформации верхней части литосферы во время рифтинга и о (главным образом) тепловом характере последующего погружения бассейна. На значения 6int при моделировании не накладывались никакие условия, т.е. полученная оценка является независимой.

Отношение интегральных факторов ßint для профилей SL и LK, равное 1.33 (точнее, его заметное отличие от отношения 5int), свидетельствует о разнице во внутреннем строеним литосферы между этими двумя сечениями. Профиль LK проходит в пределах Драбовско-Ходмской зоны фундамента с архейско-раннелротерозойским возрастом складчатости. Профиль SL располагается в районе Криворожско-Кременчугской субмеридиональной зоны интенсивной ранне-среднепротерозойской складчатости. Очевидно, что эта зональность оказывает влияние на упругие характеристики литосферы, так как древняя зональность сказывается на форме краевых разломных зон ДЦВ. Влияние это, вероятно, не слишком велико и им можно было бьг обоснованно пренебречь при менее чувствительном способе моделирования. В целом соотношение факторов растя-

жения свидетельствует в пользу гипотезы "шарнирного" характера раскрытия ДДВ (предложенной A.B. Чекуновым).

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

Информативным параметром является распределение (а не абсолютные величины) напряжений по геологическим периодам. В этом аспекте полученные диаграммы имеют следующие основные черты: сжатие в среднем карбоне; переход от сжатия к растяжению на границе карбона и перми; усиление растяжения в Перми; затем убывание растяжения до неразличимых значений.

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

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

Наиболее яркой и устойчивой чертой диаграмм палеонапряжений является "пермский максимум растяжения". Он наступает раньше, чем глобальный максимум растягивающих напряжений на границе перми и триаса, и действие его прослеживается дольше. Примечательно, что период времени между рифтингом в ДЦВ и расчетным "пермским максимумом растяжения" составляет около 60 млн. лет - что соответствует константе тепловой релаксации по МакКензи. Одним из вероятных объяснений такой картины распределения палеонапряжений во времени представляется относительно недавно сформировавшаяся гипотеза перехода от пассивного к активному механизму рифтинга.

Авалю хкороапт осадкотттстш, В связи с исследованием гипотезы перехода от пассивного к активному механизму рифтинга были проанализированы диаграммы скоростей осадконакопления, рассчитанных с учетом литофикации. На диаграммах отчетливо выделяется максимум скоростей осадконакопления вблизи оси рифта. После "пермского растяжения" максимум практически исчезает (длина волны изгибно-упругого отклика литосферы стала значительно больше 100-200 км).

"Расллывание" максимума скоростей осадконакопления со временем вполне закономерно - оно является результатом остывания ли- 17 -

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

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

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

4.4 Принципы интерпретации результатов моделирования

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

В этой связи вопрос о единственности получаемого решения особенно важен. К сожалению, как уже говорилось, в данном случае нельзя обоснованно рассчитывать на единственность получаемого решения, по крайней мере по отношению к профилю фактора /3. Практическим выходом из положения может стать применение регуляризирующих условий. Использованное в настоящей работе регуляризирующее условие придает решению устойчивость, но вопрос о его единственности в строгом смысле остается открытым. Внутри же области решений, ограниченной регуляризирующим условием, устойчивая сходимость алгоритма, вероятно, достаточна для возможности его практического применения.

Имеются и другие трудности. Как уже говорилось, стратиграфическое заполнение рифтового бассейна можно рассматривать как запись реакции литосферы на инициирующее тектоническое воздействие. Можно переформулировать это утверждение следующим образом: стратиграфическое заполнение рифтового бассейна является результатом преобразования (кодирования) первоначального (тектонического) сигнала, функцией преобразования является характеристика преобразователя-

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

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

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

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

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

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

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

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

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

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

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

На первом (исследовательском) этапе разработанный способ предоставляет возможность эффективного исследования характеристик и

- 20-

поведения модели при варьировании значений одного или нескольких параметров. Задачей исследовательского этапа является построение модели с максимальной возможной детальностью и надежностью. Модель для целей прогноза нефтегазогенерационного потенциала достаточно представить в виде 4-мерного (х, Т) поля температур и набора 3-мерных (х,1, /) поверхностей погружения для каждого осадочного слоя профиля. Модель должна быть согласована со всеми имеющимися данными уже на исследовательском этапе. Располагая такой моделью, легко восстановить палеотермобарические условия для произвольного момента времени любого стратиграфического горизонта - т.е. полную РТТ-историю его формирования.

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

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

£_ПЕИЛОЖЕНИЕ

Таб. 1 Обозначения параметров, используемых в уравнениях.

Оёозва неа и£ Нашешвзше. Ведшиш

К Мощность литосферы »125 км

К Мощность корового слоя литосферы « 45 км

Ьщ Мощность мантийного слоя литосферы к 80 ш

Мощность слоя осадочного заполнения

Рт Плотность мантийного слой литосферы я 3.3В г/см3

Рс Плотность корового слоя литосферы » 2.80 г/см3

Pf Плотность слоя осадочного заполнения

к Коэффициент температуропроводности 0.78*10"' м£/с

1 Упругий отклик литосферы на нагрузки

О(х) Изгибнэя жесткость литосферы

N Горизонтальные нагрузки, действующие на литосферу

v Коэффициент Пуассона 0.25

Е Модуль Юнга Ра

Г Температура

? Время

X Горизонтальная координата

г Вертикальная координата

/ Номер итерации

У Номер блока в модели

5го и/ Полное погружение (сумма всех мощностей осадков)

Пострифтовое погружение (сумма всех мощностей пострифтовых осадков)

Таб. 2 . Обозначения и величины параметров блоковой модели.

Оба зяяче те. Ышшаошнт Велшина

Синтетич. модель Реальные модели

1К 51

8(х) Фактор растяжения корового слоя литосферы (в каждом блоке) тах 1.98 тах 1.8 тах 1.8

Р(х) Фактор растяжения мантийного слоя литосферы (в каждом блоке) тах 2.5 тах 5.0 гпах 7.0

^ппд Время начала рифтинга (до настоящего времени) 386 млн. лет 386 млн, лет 386 млн. лет

д^пГипд Длительность рифтинга 36 млн. лет 36 млн. лет 36 млн. лет

ь Возрасты стратиграфических границ, имеющихся в разрезе см. текст см. текст см. текст

ТСЕТ Изотерма определения ЕЕТ 220 "С 220 "С 220 "С

песк Глубина некинга литосферы 10 км 10 км 17 км

N. Горизонтальные напряжения, действующие на литосферу в к-ое время см. текст см. текст см. текст

Л/ь Количество блоков в модели 200 200 200

Дх Ширина блоков 2 км 2,5 км 2.5 км

Та Температура кровли астеносферы 1333 °С 1333 ° С 1333 и С

Ьс Мощность корового слоя литосферы 40 км 40 км 40 км

К Мощность литосферы 175 км 175 км 175 км

. Рс Плотность корового слоя литосферы 2.80 г/см-5 2.80 г/см3 2.80 г/см5

Рт Плотность мантийного слоя литосферы 3.33 г/см3 3.33 г/см-" 3.33 г/см-3

а Коэффициент теплового расширения 3.4*10° "С"' 3.4*10° ° С" 1 3.4-10"3 0 С' 1

К Коэффициент температуропроводности 0.78*10"' 0.78*10'' и1/с 0.78*10"' и1/с

Л Скелетная плотность осадков 2.70 г/см3 2.70 г/см3 2.70 г/см3

Фо Исходная (поверхностная) пористость осадков 0.55 0.55 0.55

с^ Константа экспоненциального уплотнения осадков с глубиной 0.55 км"' 0.55 км"' 0.55 км '

СПИСОК ПУБЛИКАЦИЙ:

1) Поплавский К.Н. Результаты регионального автоматизированного сейсмогравитационного моделирования при изучении ДДВ. В: Проблемы повышения эффективности геологоразведочных работ на нефть и газ в XII пятилетке: Науч.-техн. конф. молодых ученых и спец.: тез. докл. - Чернигов, Б.и., 1986, - 248 е., с. 175-177.

2) Поплавский К.Н. Результаты профильного автоматизированного сейсмогравитационного моделирования при изучении ДДВ. Материалы I республиканской школы-семинара молодых геофизиков Украины, г.Алушта: Ин-т геофизики АН УССР, Киев, 1987, - 201 е., с. 81-82.:-Деп. в ВИНИТИ 5.11.87, N 7768-В87.

3) Гурова Й.Ю,, Поплавский К.Н., Козленко В.Г. Новые данные о строении скоростных разрезов Днепровского грабена. Нефтяная и газовая промышленность,: Киев, 1987, N 1, с. 9-12.

4) Козленко В.Г., Костюкевич A.C., Койфман Л.И., Кореневич К.А., Поплавский К.Н. Осадочные бассейны с высокоскоростными включениями в разрезе. В: В.И.Старостенко и Я.Шванцара (редакторы), Сейсмогравитационное моделирование при исследовании литосферы. Киев, Наук. думка,1994, с.211-240.

5) Баранова Е.П., Гурова И.Ю., Поплавский К.Н. Основные особенности скоростно-плотностных моделей ДДВ и их геологическая интерпретация. В: В.И.Старостенко и Я.Шванцара (редакторы), Сейсмогравитационное моделирование при исследовании литосферы. Киев, Наук, думка,1994, с.178-201.

6) Starostenko, V.l., Danilenko, V.A., Poplavskii, K.N., and Vengrovitch, D.B.. Modelling the origin and evolution of rift sedimentaiy basins: prc-rift subsidence. TERRA abstracts, abstract supplement N 1 to TERRA nova, v. 7, 1995, p. 61.

7) Старостенко В.И., Даниленко B.A., Венгрович Д.Б., Поплавский К.Н. Моделирование эволюции осадочных бассейнов с учетом структуры природной среды и процессов самоорганизации. Доповцц Академп наук Украши, Киев, 1996, N2, с. 97-101.

Поплавсъкий К. И. "Зворотне заедания моделювання тектокичноТ евалюци ргфтових басейтв". Дисертацш на здобуття вченого ступеня кандидата ф1зик0-математичних наук по спещальносп 01.04.12 - геоф1зика. 1нститут геофпики ш. С. I. Субботша НАН Украши, Кшв, 1996.

Захищаеться 3aci6 рйиення зворотного завдання чюельного дво.шрного моделювання процеав р^фтшга на основ1 анал1зу страттграф1-чного розр1зу, придатний як для дослщницького, так i для практичного моделювання. 3aci6 надае штерпретатору ефективний шструмент дослвд-ження модельних уявлень. Використання запронованого засобу дозволяе штерпретатору зберегги контроль за реал!зациям1 моделей, що поширюе потенцшт можливосп моделювання при доданш нових napaMeTpiB та процест.

Застосування засобу до реальних моделей Дншровсько-Донецьког западини навело до одержання ново! шформаци, що посвздчуе на користь гшотези переходу вад пасивного до активного р^фтшгу, i яка може служите додатковими пщтвердженнялш гшотези про "шаршрний" характер розкриття ДДз. Застосування засобу щодо практичного моделювання тек-тошчно! icTopii басейну дозволяе вдосконалгги прогноз нафтогазогенера-цшного потеншалу пласта.

Клюлевтлжша; модель, алгоритм, тектоничне моделювання, чисельне моделювання, р1фтшг.

Poplavskii K.N. "Inverse problem of the modeling of tectonic evolution of rift basins". Thesis for a degree of Candidate of science of Physics and Mathematics with speciality 01.04.12 - Geophysics. S.I.Subbotin Institute of Geophysics of National Academy of Sciences of the Ukraine, Kiev, 1996.

A new method for automatic searching of the best fitting parameter set of a 2D basin formation model is presented. Weighted goal functions are used in the minimization process by the inverse problem solver algorithms. The method was tested on synthetic data with parameters that were perturbed beforehand. The algorithms are capable of resolving such perturbations although there are some restrictions. The inversion technique provides the effective tool for the investigations of the rifting model's constrains. The inversion technique opens the way for further improving of the forward modeling algorithms by means of including of processes previously ignored.

The inversion technique has been applied to two profiles crossing the Dniepr-Donets Basin (Ukraine). The basic features of stress variation through time help to explain the enlarged thicknesses of Carboniferous sediments typical for this basin. The stress patterns in general agree with those derived from pa-laeo-tectonics reconstructions at both local and regional scales. However, the onset of the tension maximum in the Permian slightly precedes the global scale extensional event at the Permo-Triassic boundary and lasts longer. This may be caused by a transition from the dominance of a passive rift mechanism to an active one. The difference between two profiles modeled proves the scissors-like style of the DDB opening.