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

Российская Академия Наук Сибирское отделение Вычислительный Центр

На правах рукописи

Воронина Татьяна Александровна

• Оптимизационный подход в решении обратной задачи вертикального сейсмического профилирования

05.13.16 - Применение вычислительной техники, математического моделирования и математических методов в научных исследованиях

Автореферат диссертации на соискание \'ченой степени

кандидата физико-математических наук Новосибирск-1996

Работа выполнена в Вычислительном центре СО РАН.

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

кандидат физико-математических наук В.А.Чеверда Официальные оппоненты:

доктор физико-математических наук А.Л.Бухгейм кандидат физико-математических наук В.Г.Хайдуков

Ведущая организация: Красноярский ВЦ СО РАН

Защита диссертации состоится "23" апреля 1996 г. в " 1430" часов на заседании специализированного совета Д.002.10.02 в Вычислительном Центре СО РАН по адресу: 630090,Н°восибирск, пр. академика Лаврентьева 6 , ВЦ СО РАН

С диссертацией можно ознакомиться в библиотеке Вычислительного Центра СО РАН. Автореферат разослан " ЯЛ " -¿^ри^Ц.996 г.

Ученый секретарь специализированного совета, кандидат технических

наук ^ Г.И.Забиняко

Общая характеристика работы

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

При этом одним из наиболее развиваемых в настоящее время направлений является применение методов оптимизации (отыскание точки минимума целевого функционала) для решения обратных задач геофизики. Это обусловлено возможностью на каждом этапе численных расчетов учесть априорную информацию о строении среды и, тем самым, повысить устойчивость и достоверность получаемых результатов. Первые же попытки применения оптимизационного подхода Bamberger A., Chavent G., Lailly Р.(1979) для сейсмических приложений выявили основные трудности метода: проблема выбора начального приближения, вопросы существования, единственности и устойчивости точки минимума целевого функционала, необходимость многократного решения прямой задачи. Последнее требует наличия быстродействующих алгоритмов расчета геофизических полей.

Цель работы. Реализация подхода к решению обрати W V

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

Научная новизна работы состоит в следующем:

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

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

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

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

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

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

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

Апробация работы. По мере получения результаты диссертационной работы докладывались на Всесоюзной конференции "Математическое моделирование в геофизике" (г. Новосибирск 1988 г. ), на научных семинарах Вычислительного центра СО РАН под руководством академика А. С. Алексеева, на семинаре Института геологии и геофизики СО РАН (1993 г. ), на VI Международном Симпозиуме по Компьютерной Томографии (г. Новосибирск, 1993).

Структура и объем диссертации. Диссертация состоит из введения, трех глав, заключения и списка литературы из 70 наименований. Общий объем работы составляет 101 страницу машинописного текста, включая 23 страницы с рисунками.

Пу б л икации. О сновные результаты диссертации опубликованы в 5 работах, список которых помещен в конце автореферата.

Краткое содержание работы

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

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

ной функции и исследуется поведение целевого функционала.

Пусть полупространство {z > 0} заполнено вертикально неоднородной средой со скоростью распространения волн c(z), c(z) £ С2(0, оо), причем c(z) = с^ при z > Н. На это полупространство нормально падает плоская волна с импульсом f(t) (f(t) = 0 при t < 0 и t > Т). Волновой процесс в полупространстве z > 0 в этом случае описывается следующей краевой задачей (после выполнения преобразования Фурье по времени t):

d U + w2n2(z)u = 0, (1)

dz2

du dz

= .Р(ц;); и = Р(ш)3(ш) ехр(шс^г), г > Я, (2)

г=0

где "медленность" тг2(г) = с~2(г), а Р(ш) есть преобразование Фурье функции /(¿).

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

= и(Х}-,ш)] ] = 1,- • и>1 < и> < ш2, (3)

где (Шх,Ш2) — диапазон временных частот, на котором происходит возбуждение и регистрация волновых полей (как правило, ~ 10 • 27г, ш2 ~ 100 • 27т).

Функции п2(г) предполагаются принадлежащими множеству

М = {п2(г)еС2(0, оо); п2(г) = п^, г > Н}, на котором вводится скалярное произведете

оо

(:п\,п22)м = !(п\(г) - п22оо)(п22{г) - п^) ¿г,

после чего замыкание М становится гильбертовым пространством.

Гильбертово пространство гс2);^ 6 Ь2(и>1,и>2)',

п2 € М} будем называть пространством моделей.

Оператор п2](ш), переводящий элемент п2) этого пространства в решение краевой задачи (1)-(2) в точке Zj , назовем оператором решения прямой задачи. Очевидно, что

В, : Ь2х М —> Ь2{ш1,ш2),

то есть действует из гильбертова пространства в гильбертово.

Таким образом, если функция и^,и>) - решение задачи (1)-(2), то

£Л^;п2](и/) =

Решение обратной задачи (1)-(3) ищется как точка минимума функционала (штрафной функции)

ГШ 2

т«ч = Е / к-м - "Чиг*". (4)

3=1

В пункте 1.2 исследуется вопрос о дифференцируемо-сти по Фреше целевого функционала (4). Сначала доказывается существование сильной производной оператора В$ и выводятся формулы для нее. Достаточно показать существование частных производных по Фреше ^-[¿».^(и;) и

Выражения для градиента целевого функционала тогда имеют вид:

ль

(V* Ф^, гг2]) (ы) = -2Ке £ [«,-(«) - Е{ш)0(гЛ 0; и)] <?(*,-; ш),

¿=1

^2Ф[^,п2])(О;) = -21т [и» -

3=1

¿=11

х [«¿(о;) - Р(ш)С0; а;)] <?(*; 0; ш)<3(г-, и)йи

(5)

Здесь за (7. (г, обозначена функция Грина

задачи (1)-(2).

В пункте 1.3 исследуются стационарные точки целевого функционала. Очевидно, что он достигает минимума, равного нулю, на паре {.Р*,п2}, соответствующей истинному импульсу в падающей волне и истинному скоростному строению среды. Если он не имеет других стационарных точек, то естественно ожидать, что, по крайней мере, для достаточно "хорошего" начального приближения итерационный процесс его минимизации будет сходиться к {Г*,й1}. Пусть существует точка п|) € Ь2{ш\,и>2) X М, в которой градиент функционала обращается в ноль. Тогда из (5) следует, что

п = Еу^о^М, (6)

В диссертации доказывается

Теорема. Пусть полупространство г > 0 заполнено вертикально-неоднородной акустической средой со скоростью распространения волн с2(г:), являющейся дважды непрерывно дифференцируемой функцией, причем,

при г > Н с*(г) = Соо. Пусть в этой среде возбуждается волновой процесс, описываемый краевой задачей (1)-(2) с функцией 6 .¿/2(^1, о>2) > зарегистрированный во внутренних точках среды 21,..., гщ, Тогда градиент штрафной функции Ф[-Р; п2](п2(г) = с~2(г)), задаваемый соотношениями (5), обращается в ноль только в точке{Р8(и), такой, что

; = 1,...^, Ш е (Ш1,Ш2).

Следствие. Если в условиях предыдущей теоремы скоростное строение среды известно на некотором интервале (0, к), к < Н, причем и. по крайней мере два из цопадают в интервал (0, К), то штрафная функция ) имеет только одну стационарную точку совпадающую с истинной формой входного импульса и "медленностью" распространения волн в среде пЦг)}.

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

и

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

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

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

при ^-+1, 3 - 1,..., N - 1,

при г >

а функция и (¿; ш) удовлетворяет следующей граничной задаче

; = 1,..., ./V - 1

М

I. J г=г

Па

ди ~дг

0; 3 = 2,..., ./V, 0;

и(г,ш) = А(и>) ■ ехр(го¡гп^), при ^ > ггдт, (10)

ди ~дг

= Пы). (11)

г=0

Здесь и далее обозначено п,- = —.

Таким образом, прямая задача заключается в вычислении функции и(г, ш) по известным величинам п^,] — 1,..., N и функции

Под обратной задачей понимается определение вектора п = {п1, п2,..., п^} по известным функциям и0«(о;) = и(^,и>),б = 1,2,... ,N3. В одном варианте обратной задачи функция Р(и>) будет полагаться известной, в другом варианте будет также подлежать определению.

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

Для вычисления градиента штрафной функции (4), как следует из (5), необходимо знание функции Грина 0(г,£,ш). Под функцией Грина будем понимать решение задачи (8)-(11) при ^(ш) = 1. Тогда решение задачи (8)—(11) будет получаться умножением функции Грина на F(u). В силу принципа взаимности функцию Грина О (г, со) для этой модели можно найти как решение краевой задачи, отличающейся от задачи (8)—(11) только тем, что условие (11) заменяется

условиями:

Эи Эг

ди дг

г=0

= 1.

(12)

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

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

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

Целевой функционал, подлежащий минимизации, имеет вид

м, у

= £ / (13)

5=1 £

Здесь имеется в виду, что присутствует зависимость подынтегральной функции от вектора {п2, к = 1,...,ЛГ}, а именно, от этих значений зависит функция Грина 0(г:,ш) где 0(гВ) 0, си) = С(23)ш) — 2-ой аргумент для краткости опускаем.

Градиент функционала п2) тогда является вектором,

к- компонента которого имеет вид:

дФ ' ^

дп\

= 2Ке /£{(Р(си)С(28,а;)-и8(а;))^^Р(а,)}аа;. 8=1 к

(14)

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

• при каждом фиксированном ш по алгоритму, описанному в предыдущем разделе, вычисляется набор значений {(7(^,0;), й = 1,.. . и

• после чего вычисляются значения подынтегральных функций в формулах (13)—(14);

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

Далее рассматривается задача с неизвестным импульсом. Теперь определению подлежат вектор {тг|, к = 1,..., iV} и спектр импульса F{ui){u 6 [u>i,u2]}, доставляющие минимум целевому функционалу (13).

Поскольку этот функционал зависит от F(u>) квадратичным образом, то минимизация по F(w) может быть проведена аналитически и, тем самым, мы исключим из задачи минимизации F(u>) и преобразуем целевой функционал к но-

ФЫ = 1 {"^мГ+S кН| )^

(15)

В этом выражении легко видеть величины, участвующие в классическом неравенстве Коши-Шварца. Имеется квадрат скалярного произведения двух ^-мерных векторов и квадраты их модулей. В силу этого неравенства функционал равен нулю только тогда, когда при каждом ш векторы {us, 5 = 1,..., Ns} и {Gs, 5 = 1,..., Ns} пропорциональны друг другу. То есть, при каждом ш существует такой коэффициент F, что иs = F ■ Gs для всех s. Этот коэффициент и является значением спектра импульса F(uj).

Таким образом, задача оптимизации теперь решается без обращения к функции F(u>), причем при желании на любом шаге можно узнать расчетный импульс F(us) по формуле :

^И - г Л

Глава 3. Третья глава состоит из трех пунктов, в которых проведен подробный анализ численных результатов.

В п о. 1 обсуждаются вопросы выбора начального приближения среды. Как известно, одной из основных проблем здесь - реконструкция трендовой составляющей, т. е. плавной "огибающей" скоростного разреза. Впервые эта проблема возникла еще в работах Kunetz G.(1961,1963). При этом оказалось, что восстановление идет достаточно устойчиво, если в спектре входного сигнала присутствуют временные частоты, сколь угодно близкие к нулю. При реализации оптимизационного подхода возникла схожая проблема:

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

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

В п 3. 2 рассматриваются оптимизационные методы поиска точки минимума целевого функционала. Минимизация осуществляется методом сопряженных градиентов, а поиск

минимума на каждом шаге осуществляется методом золотого сечения.

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

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

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

На приведенном рис.1 представлен характерный результат восстановления среды при неизвестном импульсе. Мож-

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

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

и. У

СЬ)

аз

аб

1.0

ч.

№

г

0.0

0.3

0.Б

[(7)

0.0 0.7 1.4 2.1 0.0 0.7 1.4 2.1

а) 0.«

б)

в)

0.0 0.7 1.4 2.1

Рис. 4 Результат восстановления скоростного строения среды. Дополнительный приемник на глубине 140 м (диапазон временных частот от 10 до 160 Гц): а) - истинная скорость; б) - начальное приближение; в) - результат восстановления.

В Заключении сформулированы основные результаты диссертации.

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

Работы автора по теме диссертации.

[1] Воронина Т.А., Чеверда В.А. Модельная обратная задача вертикального сейсмического профилирования.-Новосибирск, 1993.-52с.-(Препринт / РАН. Сиб. отд. ВЦ;1004).

[2] Воронина Т.А., Чеверда В.А. Обращение полных волновых полей при обработке данных метода сейсмического профилирования // Доклады Академии Наук.- 1994.-Т. 335, №4.

[3] Воронина Т.А., Чеверда В.А. Оптимизационный подход к обработке данных метода вертикального сейсмического профилирования // Геология и геофизика-1994. №5, С.127-139.

[4] Cheverda V.A., Voronina Т.А. Full waveform inversion for seismic borehole date //J. of Inverse and Ill-Posted Problems, 1994.-V.2, №3,-pp.211-226.

[5] Cheverda V.A., Voronina T.A. Full waveform inversion of VSP data (normally incident plane wave). Bull.Nov.Сотр.Center: Math. Modelling in Geophysics, 1994, n.l.