О новом подходе к решению обратной динамической задачи сейсмики для слоистых сред

О новом подходе к решению обратной динамической задачи сейсмики для слоистых сред

А.С. Рожков

НИИ Моргеофизика

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

Abstract. The existence and the uniqueness of the solution of the inverse dynamic seismic problem have been proved and the way of its solution for the layered medium under rather common assumptions about its features has been given in the paper.

1. Введение

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

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

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

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

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

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

Опишем формально требования, налагаемые на среду. Отложим по горизонтали ось ОХ с началом в некоторой точке О. Положительным направлением будем считать направление слева направо. Пусть имеется п слоев, для каждого из которых в некотором отрезке [х0, X] значений х известны жесткости С,(х) и времена пробега сейсмических волн /¿(х), как функции от х, т.е. хе[х0, X], i - номер слоя, i = 1,...,n. Число границ, которые разделяют отдельные слои, равно n+1. Пусть известны глубины первой и последней из этих границ как функции от х, т.е. И1(х) и Н„+1(х). Обозначим мощность ¿-го слоя через Н,(х). Тогда

Z Н

Поставим перед собой задачу найти в каждом i-м слое функции Н(х) - мощность, У(х) -скорость, Р(х) - плотность слоя. Зная И1(х) и найдя Н(х), тем самым найдем и неизвестные И(х),

/ = 2,...,п. Мы предполагаем, что У(х), Р(х) являются для всякого фиксированного значения хе[х0, X] постоянными по глубине в пределах / -го слоя.

3. Первый способ решения задачи

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

С

где / - номер /-го слоя, / = 2,...,п.

Кроме того, из работы (Гогоненков, 1987) следует, что удовлетворительная аппроксимация зависимости Р от У может быть получена при использовании уравнения

Р = тУ п,

где т для разных районов изменяется от 0,293 до 0,902, а п - от 0,134 до 0,902. Из уравнения (3) получим Р1/п = т1/пУ, поэтому У = т'1/пР1/п. Положим т'1/п = с, 1/п = I, тогда

У = cPl,

3,73 < I < 7,46.

В работе (Гогоненков, 1987) говорится также, что для основных литотипов пород (песчанники, глины, известняки, доломиты) связь между У и Р удовлетворительно описывается уравнением (4). Гидрохимические осадки (соль, гипс, ангидрит), а также уголь характеризуются аномальным соотношением У и Р. Для некоторых пород, например, известняков, зависимость близка к линейной. Из уравнений (2) и (4) следует, что

С = сР1+1.

С'=

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

Производная этой функции непрерывна всюду на отрезке [х0, X], кроме тех точек, которые являются точками соединения вышеупомянутых отрезков. В силу уравнения (6') производная функции С(х) непрерывна в тех же точках, что и функция Р(х), а в тех точках, где производная функции С(х) разрывна, терпит разрыв и функция Р(х). Зная функцию С(х), мы можем определить точки, в которых ее производная терпит разрыв. Тем самым мы найдем и все точки, в которых соединяются отдельные звенья функции Р(х).

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

Итак, пусть мы разбили отрезок [х0, X] на отрезки точками х-,- = 0,...,к такими, что выполняются

условия

х0<х!<.<хы<X, хк = X.

Возьмем некоторый из отрезков [х-х-н],- = 0,...,к-1. Введем систему координат с началом в х-. На каждом из этих отрезков Р(х) является просто линейной функцией вида

Р

где Р- - значение плотности в начале --го отрезка, а- - некоторые константы, т.е. Р = Р(х-), - = 0,...,к-1. Если перейти теперь в систему координат с началом в точке х0, то для --го отрезка Р(х) = а-(х-х-)+Р-, хе [х-, х-+1]. Из (4) и (6) будем иметь

У

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

плотностью соответственно уравнениями (4) и (6), то скорость и плотность имеют производные первого и второго порядков в тех же точках, что и плотность.

Из (9) получим:

У

У"

Из (9), (10) и (11) получим:

1 = У(х)У"(х)/((У(х))2 - У(х)У'(х)) + 1,

I = С

[а)(х-х)) + Р)]/а) = У(х)У(х) / ((У(х))2- У(х)У"(х)),

[а)

из (13) будем иметь:

из(14) получим:

Р) /а) = У(х)У(х)/((У(х))2- У(х)У'(х)) - (х-хД

Р)/а) = С

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

Мы можем знать значения функций У(х) и С(х) лишь в дискретном наборе точек. Поэтому, чтобы вычислить У(х) и У"(х) или С(х) и С"(х) необходимо применить численные методы. Эти методы широко известны. В работах (Березин, Жидков, 1962; Бахвалов и др., 1987; Самарский, Гулин, 1989) изложена подробная теория и даны практические алгоритмы нахождения производных. Даны также оценки вычислительной погрешности при вычислении производных. В данной работе мы не будем останавливаться на этих вопросах.

Итак, на любом)-м отрезке ) = 0,...,&-1 по формулам (15) можем найти Р)/а), где а) = Р'(х) на)-м отрезке, а Р) - значение Р(х) в начале отрезка.

Заметим, что выражения (15) справедливы для всех точек х из отрезка [х), х^]. Из этого следует, что если мы найдем по формулам (15) (приближенно) значения Р)/а) для нескольких точек отрезка [х), х)+1], то затем можно вычислить среднее значение этих, приближенно вычисленных, величин. Так как истинное значение Р)/а) единственно (в силу наших предположений о линейности функции Р(х) на каждом отрезке) для всего этого отрезка, то, как следует из теории вероятностей, при увеличении числа точек х, для которых мы по формуле (15) вычисляем отношение Р)/а) (с некоторой погрешностью), среднее значение этих величин стремится к истинному. Следовательно можно значительно уточнить отношение Р)/а) для каждого отрезка.

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

Итак, имеем:

Р(х) = Р0 + а0(х-х0) на отрезке [х0, хх],

Р(х) = Р\ + а!(х-х!) на отрезке [х1, х2],

Р

Из (16) получим:

Р(х) = Р)[1 + (а)/Р))(х-х))] на отрезке [х^, х)+1]. (17) Подставляя в (17) вместо х его конечное значение х)+1 на отрезке [х^, х)+1], и, имея в виду, что Р(х)) = Р) получим:

Р\ = Р0[1 + (ао/Р0)(х!-х0)] на отрезке [х0, хх],

Р2 = Р1 [1 + (а\/Р\)(х2-х1)] на отрезке [хь х2],

Р3 = Р2[1 + (а2/Р2)(х3-х2)] на отрезке [х2, х3],

Р) = Р)-1 [1 +

Подставляя во второе из уравнений (18) вместо Р\ его значение, получаемое из первого уравнения, получим:

Р2 = Р0[1 + (а0/Р0)(х!-х0)][1 + (а^РО^-хО] на отрезке [хь х2]. Теперь в третье из уравнений (18) подставим вместо Р2 его значение, стоящее в правой части последнего уравнения. Получим:

Рэ = Р0[1 + (а0/Р0)(х:-х0)][1 + (а:/Р:)(х2-х1)][1 + Ц^^^)] на отрезке [х2, хэ]. Поступая аналогично с последующими уравнениями (18), придем к рекурентной формуле:

Р- = Р0[1 + (а0/Р0)(х:-х0)][1 + (а]/Р])(х2-х|)]...[1 + (ая/Ря)(х--хя)] на отрезке [хя, х-]. (19) Для удобства введем теперь следующие обозначения:

Е0= 1, Е- = [1 + (ас/Р0)(х1-х0)][1 + (а^РО^-хО]...^ + (а]Л/Р]Л)(хгх]Л)], -=1,...,к-2. (20) Тогда из (17), (19) и (20) получим

Р

Введем теперь функцию/х), такую, что для всех хе [х-, х-+1]

Ах) = Е-[1 +

Эта функция определена во всей области определения функции Р(х) и, т.к. числа Е-, ц/Р-, ху -известны, то Ах) - известная функция на всем отрезке [х0, X]. Из (21) и (22), наконец, следует, что

Р

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

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

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

Сначала найдем методом, изложенным выше, в каждом слое функции Р,(х) с точностью до постоянного множителя Р,(х0), т.е. по формуле (23) получим:

Р,

Из (2) следует, что

У

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

Н

а из (24) будем иметь

Н

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

Л, = 1/Р,

Заметим, что О,(х), в силу сказанного выше, - известные функции, а Л, - неизвестные числа для всех /=1,2,...,п. В связи с этим (27) можно записать следующим образом:

Н,

Подставив выражения, полученные в (29) для Н,(х), в формулу (1) получим:

ЕН,

Разобьем теперь отрезок [х0, X] на оси ОХ некоторым образом на п-1 частей точками хьх2,...,хп

(рис. 1).

Подставив в (30) вместо х последовательно каждое из х-,-=1,...,п, получим:

А1в1(х1)+...+Апвп(х1) = Я(хО,

А101

Получим систему линейных уравнений п-го порядка с неизвестными А,-. Найдя А, найдем из равенств (28) и Р,(х0). Затем, так как функции /(х) были вычислены ранее, найдем по формуле (24) функции Р(х). Затем по формуле (25) найдем функции У(х), а по формуле (26) - функции Н,(х). Из единственности решения линейной системы уравнений следует единственность найденных функций Р,(х), У(х) и Н,(х). Тем самым мы полностью решили сформулированную нами задачу.

4. Второй способ решения задачи

Теперь приведем другое решение той же задачи. Заметим, что при решении мы ни разу не воспользовались формулами (12). По формулам (12) можно вычислить число 1, являющееся степенью при Р в уравнении (4). Это число является одним и тем же для всего слоя. Поэтому оно может быть найдено с использованием всего множества отрезков [х^, х^ц] )=0,1,...,&, в то время, как отношение а)/Р) находится для каждого ) только по одному отрезку [х), х^]. Следовательно, точность вычисленного

значения I по формулам (12) может превышать точность вычисления чисел а/Р, по крайней мере для

1 1

небольших значений I. Пусть мы нашли число I. Из формулы (4) получим: У(х) = сР(х) и У(х0) = сР(х0). Из последних выражений получим:

У(х)/У(х0) = [Р(х)/Р(х0)]г,

откуда

Р(х)/Р(х0) = [У(х)/У(х0)]-1,

поэтому

Р

Аналогично для функций Р(х) и С(х), используя формулу (6), получим:

Р

Формулы (33) и (34) аналогичны формулам (23). Но, в отличие от них, эти формулы могут оказаться значительно точнее, так как число I находится точнее, чем числа а)-/Р)-. К тому же для вычисления по формулам (23) (если I не велико) производится значительно больше операций, чем при вычислениях по формулам (32). А ошибки вычислений при использовании величин имеющих погрешности, возрастают с ростом числа операций, произведенных над этими величинами.

С другой стороны, первое решение обладает тем преимуществом, что не требует вообще знания числа I. Это означает, что если в отдельном слое число I не остается постоянным, то это не влияет на результат вычисления.

Таким образом, для каждого слоя мы можем найти по формулам (33) и (34) функции Р,(х) с точностью до постоянного множителя в виде Р,(х) = Р,(х0)/(х).

Повторяя далее ход решения первого решения задачи, который следует за выводом формул (23), мы окончательно решим задачу.

5. Заключение

Скажем здесь несколько слов по поводу возможности решения линейной системы (32). Теория решения линейных алгебраических уравнений является хорошо разработанной теорией. Полно и всесторонне эта теория освещена во множестве работ разных авторов. Одними из наиболее значимых работ по освещению этой теории, по мнению автора, являются работы (Воеводин, 1977; Хорн, Джонсон, 1989; Тихонов, Арсенин, 1986).

В работе (Хорн, Джонсон, 1989) на стр. 406 приведена оценка относительной ошибки в решении линейной системы типа (32), содержащей ошибки в начальных данных. Не приводя здесь этой оценки, приведем высказывание из этой работы на той же странице. "Число обусловленности %(А) ... играет решающую роль и определяет чувствительность оценки ошибки в решении к погрешностям в начальных данных".

Число обусловленности матрицы А по отношению к матричной норме ||А|| определяется следующим образом:

если А не вырождена,

%(А) = »,

если А вырождена.

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

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

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

Литература

Алексеев A.C. Обратные динамические задачи сейсмики. В кн.: Некоторые методы и алгоритмы

интерпретации геофизических данных. М., Наука, с.9-83, 1967. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М, Наука, 600 с., 1987. Березин И.С., Жидков Н.П. Методы вычислений. М, Государственное издательство физико-

математической литературы, т.1, 464 е., 1962. Воеводин В.В. Вычислительные основы линейной алгебры. М., Наука, 188 е., 1977.

Гогоненков Г.Н. Изучение детального строения осадочных толщи сейсморазведкой. М., Недра, с.221, 1987.

Рожков A.C. Решение обратной динамической задачи для слоистых сред. Вестник МГТУ, т.3, № 1, с.79-80, 2000.

Рожков A.C. Способ автоматического определения точек излома вертикального годографа ВСП.

ВИНИТИ РАН, № 2268-В99, с.8, 1999. Самарский A.A., Гулин A.B. Численные методы. М., Наука, 432 е., 1989. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М., Наука, 288 е., 1986. Хорн Р., Джонсон Ч. Матричный анализ. М., Мир, 655 е., 1989.