CC BY
58
Определение трехмерной слоисто-неоднородной скоростной модели среды по данным метода отраженных волн
Э.А. Бляс
Судоводительский факультет МГТУ, кафедра высшей математики
Аннотация. Рассматривается задача определения трехмерных скоростных моделей по данным сейсморазведки методом отраженных волн. Дается описание основных методов решения данной задачи, делается сравнительный анализ, рассматриваются преимущества и недостатки каждого метода. Приводится новый итерационный алгоритм определения модели среды с неоднородными криволинейными слоями, не требующий пересчета кривизн фронта волны. Показано, что для определения скоростей неоднородных слоев достаточно уметь решать обратную задачу в рамках модели с локально-однородными слоями. Приводится оптимизационный алгоритм, не требующий многократных решений прямых задач, что позволяет существенно сократить время счета.
Abstract. The given paper deals with the problem of 3-D velocity layered model definition from reflection seismic data. It gives the description of basic methods for the named problem, gives the compare analysis, discusses advantages and disadvantages of these methods. It demonstrates a new iteration algorithm for the determination of the medium model with inhomogeneous curve-linear layers which does not require the wave front curvature. It is evident from the work that in order to define velocities in the inhomogeneous strata one has to solve an inverse problem for the model with locally-homogeneous layers. It demonstrates the optimization algorithm that does not require multiple solution of the direct problem and reduces time-consumption significantly.
1. Введение
Рассмотрим задачу определения трехмерной скоростной модели среды с криволинейными горизонтально-неоднородными слоями по кинематическим характеристикам отраженных волн. Для ее решения можно применять оптимизационные методы, широко используемые при решении обратных задач геофизики (Голъцман, 1971; Гольдин, 1979; Гольдин и др., 1993; Яновская, Порохова, 1983). В работах С.В.Гольдина подробно рассмотрены различные методы оптимизации, применяемые при решении обратных двумерных кинематических задач. Они основаны на многократном (сотни раз) решении прямой кинематической задачи и минимизации расхождения между наблюдаемыми и рассчитанными временными полями отраженных волн. В принципе эти методы позволяют учесть кривизну лучей в неоднородных слоях, то есть найти скоростную модель в классе моделей с неоднородными слоями, но для их применения требуется хорошее начальное приближение. Кроме того, при наличии сильных горизонтальных неоднородностей решение прямой задачи в среде с криволинейными неоднородными слоями требует больших затрат машинного времени даже в двумерном случае. Для трехмерных моделей сред с горизонтально-неоднородными несогласными слоями затраты машинного времени существенно возрастают, но, что еще более важно, для описания поверхностей слоев и пластовых скоростей может понадобиться большое число базисных функций, что уменьшает устойчивость решения обратной задачи. Кроме оптимизационных алгоритмов, для восстановления скоростной модели среды по эффективным параметрам метода ОГТ можно применять итерационные инверсные алгоритмы (R-алгоритмы по терминологии С.В.Гольдина), основанные на приближенном решении обратной задачи (Глоговский, Гогоненков, 1978; Глоговский, 1985; Sattleger, 1975; Гольдин, 1986). Сходимость такого алгоритма для двумерных моделей слоисто-однородных сред рассматривалась в (Бляс, Левит, 1989), где показано, что в большинстве практически важных случаев он сходится за 2 - 4 шага. При построении R-алгоритмов определения скоростной модели слои предполагаются локально-однородными, так как только в локально-однородном слое можно применить инверсный метод определения его скорости. Напомним, что слой является локально-однородным, если в пределах отдельного луча его можно считать однородным, то есть пренебречь кривизной луча и связанным с ней изменением времени.
Таким образом, хотя итерационные R-алгоритмы быстро сходятся, используемые в них инверсные методы разработаны только для среды с локально-однородными слоями. С другой стороны, оптимизационные методы позволяют восстанавливать скорости с учетом локальной неоднородности слоев, но требуют больших затрат машинного времени. В настоящей работе рассматриваются алгоритмы
обоих видов. Недостатком методов оптимизации является необходимость многократного (сотни раз) решения прямой задачи, которая для реальных многократных схем наблюдений сводится к расчету десятков тысяч лучей. В средах с сильными криволинейными преломляющими границами в верхней части разреза на эффективную скорость существенное влияние оказывают вторые производные этих границ, причем это влияние увеличивается с ростом глубины залегания восстанавливаемого слоя. Это значит, что для определения скоростей глубоко залегающих слоев необходимо верхние резкие границы восстанавливать с точностью до вторых производных (Бляс, 1988, 1991). Как известно, отражающие границы по реальным данным восстанавливаются только с точностью до первых производных (если учитывать направление подхода нормального луча), поэтому с увеличением глубины погрешности построений растут. Исходя из этого, можно предложить следующий подход к определению скоростных моделей в средах с сильными горизонтальными неоднородностями в покрывающей толще. Задача решается послойно сверху вниз. Сначала применяется Я-алгоритм к определению текущей границы, а затем применяется оптимизационный метод, при котором уточняются все найденные ранее слои. Использование априорных ограничений на параметры определяемых скоростей и границ при минимизации целевой функции сразу по параметрам всех слоев позволяет уменьшить регулярные погрешности восстанавливаемых глубоких слоев за счет влияния кривизны сильных преломляющих границ в покрывающей толще, так как в таком процессе уточняются не только отражающие границы и скорости в покрывающих их пластах, но и преломляющие границы, определенные ранее. Здесь необходимо отметить, что оптимизационные методы применяются также при определении скоростной модели среды по данным метода вертикального сейсмического профилирования (ВСП), когда приемники располагаются в скважине. В этом случае понятие скорости суммирования (аналогично скоростям суммирования при поверхностных наблюдениях) неприменимо, поэтому оптимизация является практически единственным методом уточнения скоростной модели околоскважинного пространства.
2. Постановка задачи
Пусть ¥к(х,у) - функции, описывающие поверхности слоев, ук(х,у) - пластовые скорости, к -номер слоя. Предположим, что имеется несколько сейсмических профилей, покрывающих исследуемую площадь. На каждом профиле известны линии ^ - времена пробега волн с нулевым расстоянием взрыв-прием, эти времена известны для всех отражающих границ. Кроме того, предположим, что после проведения скоростного анализа известны скорости суммирования уогт. Хотя скорости суммирования определяются по сейсмограммам и зависят от волнового поля, во многих случаях можно считать, что эта скорость определяется как параметр гиперболы, аппроксимирующей годограф ОГТ (сПР(1) по методу наименьших квадратов на базе наблюдений, I - расстояние взрыв-прием (Пузырев, 1979). Скорости суммирования находятся для всех отражающих границ. Фактически количество регулярных отражений вместе со скоростями ОГТ (суммирования при получении разрезов ОГТ) и определяет количество слоев в скоростной модели. Таким образом, на каждом из профилей имеется 2п известных функций: (0к, УсаР(^к1, к=1,2,...,п. Требуется определить пластовые скорости ук(х,у) и отражающие границы ¥к(х,у).
Прежде чем рассмотреть итерационный Я-алгоритм, необходимо рассмотреть инверсный метод решения обратной задачи, применение которого является составной частью итерационного алгоритма.
3. Инверсный алгоритм определения скоростной модели
Задачу определения скоростной модели будем решать послойно сверху вниз. Сначала находим скорость и мощность первого слоя по известным кинематическим параметрам волн, отраженных от его подошвы. После этого, считая параметры первого слоя известными, находим подошву и скорость распространения волн во втором слое и т.д. Таким образом, инверсный алгоритм строится в предположении, что неизвестными являются параметры только одного слоя, а именно, того, от подошвы которого выделены отраженные волны. Обратная задача решается по эффективным кинематическим параметрам этих волн.
Итак, предположим, что известны параметры первых (п-1) слоев, требуется определить скорость у(х,у) и подошву ¥п(х,у) п-го слоя по известным линиям 4(р), у0(р), grad ?0(р) по некоторой системе профилей, р - координата точки на профиле. Для каждой точки профиля будет найдена скорость уп в окрестности центрального луча, связанного с этой точкой, и восстановлена точка отражения этого луча от подошвы слоя. Имея достаточно плотную систему профилей, можно затем сглаживанием восстановить отражающую поверхность ¥п(х,у) и скорость уп(х,у) п-го слоя. Данная задача для однородных слоев рассматривалась во многих работах (Черняк, Гриценко, 1979; Гольдин, 1993), здесь мы приводим другой алгоритм, не связанный с пересчетом кривизн фронта волны.
Пусть (Х,У) - координаты точки ОГТ на некотором профиле. Зная grad 4, можно, используя закон Бенндорфа (Гольдин, 1979), восстановить направление подхода центрального луча к точке (Х,У) и
протрассировать луч до (п-1)-ой границы - до кровли п-го слоя; эту границу будем обозначать через Р(х,у), отбрасывая индекс "п-1" для сокращения записи. Точку пересечения данного луча с (п-1)-ой границей обозначим через (^щР^, F0 = Р(^0ь^0). Далее воспользуемся равенством, впервые строго доказанным С.А.Гриценко (Черняк, Гриценко, 1979; Гольдин, 1986)
2 /(^2) = d 2. (1)
Здесь ф) - эйконал лучей, вышедших из (искомой) точки M(%r„цn,Fn) отражения центрального луча и пришедших в точку профиля, удаленную (с учетом знака) на расстояние s от точки А0(Х,У) (рис.1); ve -дифференциальная скорость, найденная по профильному годографу ОГТ с центром в точке А0(Х,У), Ь -прямая, на которой расположены источники и приемники (профиль наблюдений). Через обозначим точку пересечения с (п-1)-ой границей луча, вышедшего из точки М и пришедшего в точку профиля на расстоянии 8 от точки А0(Х,У), то есть N(5) - точка пересечения с (п-1)-ой границей луча, соответствующего эйконалу ф). Тогда ф) можно представить в виде
г= 71(5,^) + Ц^) = Г(5,^). (2)
Здесь Т1 - время пробега волны из точки до точки А на профиле, расположенной на
расстоянии 5 от А0, tn - время пробега волн п-ом слое. Так как оператор ф решения обратной задачи строится для среды с локально-однородным п-ым слоем, то в этом слое луч представляет собой прямолинейный отрезок, и для tn справедливо равенство
4 = 1Т))2+ (&-£) + (%- 7)2]2/ V. (3)
X
Из принципа Ферма (стационарности траектории) следует, что выполняются следующие равенства:
сТ /д$ = дТ^г) /д$ + дп(£,,г)) /д$ = 0,
(4)
сТ /дг] = дТ^^г) /дц + д^(£,,г) /дп = 0. Подставляя в данные равенства tn из (3), получим, что
[(£ - 4) + (Г - Fn) (РХ)] / [vn [(F - Fn)2 + а - & )2 + (Л - %)2]1/2] = - дТ /д$,
[(?- Чп) + (Г - рп) (дР/Щ / [vn [(Г - рП)2 + (£- ^)2 + (?- 7П)2]1/2] = - дТ/дт].
Пусть О = Р-Рп - неизвестная величина. Тогда из (3) и (5) получаем, что выполняется равенство
[(¿Р/^) О + (сТ^д^) уп А(]2 + [(сР/дфО + (сТ1/с^) уп А^2 + О2 = (уп А()2. (6)
Здесь и всюду ниже производные функций Т1, Р берутся соответственно в точке (0,^0), ^(0)) и (£(0),^(0)) соответствующей точке пересечения центрального луча с (п-1)-ой границей, А! = ^/2 - Т1 - известная величина. Координаты проекции точки пересечения луча, выходящего из точки М, зависят от 8, то есть являются функциями переменной 5, поэтому £(0), ^(0) означает, что функции берутся в точке 5 = 0, что соответствует центральному лучу.
Равенство (6) после несложных тождественных преобразований можно записать в виде
Б1О2 + В2 О Уп2 + БэУп4 + Б4У2 = 0, (7)
где
В1= 1 + (сР/Зс)2 + (сР/ду)2, В2 = 2(сР/Зс ЗТ1 /д$ + дР/ду сТ /дг) Мп,
Вз = (А^)2 [(¿Т1Щ )2 + (¿Т1 /дф2], В4 = - (А^)2.
В уравнении (7) неизвестными являются величины Уп и О, для их нахождения необходимо иметь еще одно уравнение. Его получим из равенства (1). Для нахождения й 2т /й5 воспользуемся способом, впервые изложенным в (Бляс, 1985), а именно, продифференцируем (2) два раза по 5 с учетом того, что Е,, ц являются функциями 5. Используя равенства (4), получим
йгй = дТ/с5, й2т/д25 = д2Т1 /¿52 + (,д2Т1 /¿5 + (д2Т1 /с5 йцЩ. (8)
Далее, дифференцируя (4) по 5, получим систему линейных уравнений относительно й£/йя, йц/йв:
{д2Т/д£,2) (й£,Ш) + (д2Т/д£,дг))(йчШ) = - дгТ1/д^с5, (д2Т/д$дг}){й$Ш) + {д2Т/д2г) (й^Щ = - д2Т1/дг1д5.
Решая данную систему относительно й^/й^', йц/й$ и подставляя их во второе равенство (8), найдем й 2т /й^2. Подставляя найденное таким образом й 2т/й52 в равенство (1), получим
2 /{0 У2) = (32Т1/с52) - {(^2Т1/^^ск) [(^¡ЩсЬ^Т/дг?-) - ^Т^дцсЪ) (д2Т/д£,дг1)] +
+{д2Тг/дг!д5) [(д2Тг/дг1д5){{д2Т/д1;2) - (с>2Т1/с>%&) (д2Т/д£,дг))] }/ [(д2Т/д£)(д2Т/дт?) - (¿>2Т/^7)2].
Подставляя в это равенство Т из (2) и вычисляя вторые производные величины tn, определенной равенством (3), получим равенство, которое можно записать в виде
А1 О2 + А2О + А3ОУ2 + А4У2 + А5Уп4 = 0. (9)
Здесь
(Ап А1 = [(32Т1/с52) - 2 /(}0 у2)] [(^2Р/Зс2)(^2Р/су2) - (д2Р/сСсу)2],
(Ап)2А2 = - [(д2Т1/с52) - 2 /(0 у2)][(д2Р/ск2)(1+Р2) + (д2Р/су2)(1+Р2) - 2(д2Р/сксУ)РхРу],
АЫз = (д2Т1/д%&)2 (д2Р/су2) + (д2Тг/дг1д5)2 {д2Р/дс2) - 2 (д2Т1/д1;д5){д2Т1/дг1д5){д2Р/дс^ду) -- [(^/Л) - 1 /(^у2)] [(^2Р/ск2)Б2 + (д2Р/су2)В1 - 2 (д2Р/сУск) Бз],
МпА4 = - (д2Т1 /д£,о5 )2 (1+Ру2) + (д2Тг/дг1д5)2(1+Р2) - 2{д2Т1/д1;д5){д2Т1/дг1д5){дР/сС){сР/ду)-- [(^2Т1 /д28) - 1 /(^у2)] [(1+Рс2)Б2 + (1+Р2)Б1 - 2 Рс Ру Бз)], АЫ5 = (д2Т1/д&)Б2+(д2Т1/дг1с5)П1 - 2(д2Т1/д&)(д2Т1/д1]с5)Оз-- [(^2Tl/йS2) - /(^ у2)](Б102- Бз2) ,
где
(10)
Di = (d2T1 2) - {дГ1 /д§2 (Atn) -1, D2 = (d2T1 /dif) - (T /árj )2 (At„ )-1, D3 = (d2T1 /dE,dr¡) - (T /^(T /drf)(Atn )-1.
Напомним, что производные функций T1, F берутся соответственно в точках (0,^0, r¡0) и соответствующих точке пересечения центрального луча с (п-1)-ой границей z = F(x,y). Итак, для нахождения неизвестных G, vn имеем два уравнения (7), (9). Пусть а = G /vn2, тогда данные уравнения можно записать в виде
B1 о2 v2 + B2 av2 + B3V2 + B4= 0,
(11)
A1 о2 V2 + A2 а + A3 av2 + A4 + A5V2 = 0.
Из первого уравнения (11) находим vn2 = - B4/(B1 а2 + B2 а + B3). Подставляя отсюда vn во второе равенство (11), получаем кубическое уравнение относительно а:
B4 (Aia2 + A3a + A5) = (A2a + A4)(B1a2 + B2a + B3). (12)
Таким образом, нахождение пластовой скорости в окрестности центрального луча сводится к решению кубического уравнения (12) и нахождению vn по формуле
vn= [- B4 /(B3+B2a+B3c¿2)] 1/2
Как аналитически показано автором (Бляс, 1991), с увеличением глубины залегания преломляющей границы влияние ее вторых производных на ve уменьшается. Это значит, что для получения приближенного значения пластовой скорости можно не учитывать влияние вторых производных кровли оцениваемого слоя, то есть при построении оператора ф полагать d2F/ck 2=d2F/3cdy =d2F/dy 2= 0. Тогда из (10) видно, что А1=А2=А3=0 и из уравнения (9) сразу находим vn:
vn= (- A4/A5)1/2
В формулы (10) входят вторые производные временного поля T1(s,^, rj) вдоль центрального луча, эти производные можно найти методом, изложенным в (Бляс, 1985). Итак, данный алгоритм позволяет найти скорость вдоль центрального луча, а значит, и точку отражения в предположении локальной однородности оцениваемого слоя. Заметим, что выше залегающие слои не предполагаются локально однородными, так как до оцениваемого слоя центральный луч и производные вдоль него можно считать с учетом кривизны лучей (Приложение 2).
4. Итерационный R-алгоритм
Рассмотрим итерационный R-алгоритм, опирающийся на изложенный инверсный метод определения пластовых скоростей и отражающих границ в трехмерных слоистых средах. Данный алгоритм также работает послойно, поэтому предположим, что известны (ранее найдены) параметры первых (n-1) слоев - функции Fk(x,y), vk(x,y), k=1,2,..., n-1. Требуется определить скорость vn(x,y) и отражающую поверхность z=Fn(x,y) n-ro слоя по известным функциям t0, vcdp волн, отраженных от этой поверхности и полученных по полевым сейсмограммам. Данную задачу можно рассматривать как задачу нахождения решения операторного уравнения
f(S)= d, d = (to,vcDp). (13)
Здесь S - искомая модель, то есть подошва Fn(x,y) и скорость vn(x,y) в слое: S = (Fn(x,y), vn(x,y)), f -оператор решения прямой задачи, действие которого сводится к следующим операциям. Сначала численными методами для заданной модели S = (Fn(x,y), vn(x,y)) слоя рассчитываются годографы ОГТ для заданного набора профилей. После этого находится интегральная скорость vCDP как параметр гиперболы, аппроксимирующей годограф ОГТ tCDP(l) по методу наименьших квадратов, то есть из условия минимума функции
I
.2 , ,2/. 2\1/2т2
ф(успр) = i [сор(и) - (0 + 1,2/Успр2)т]\ (14)
1=1
где и - расстояние взрыв-прием на /-ой трассе сейсмограммы ОГТ, I - количество точек в годографе ОГТ.
Уравнение (13) означает, что необходимо найти такие значения пластовой скорости У„(х,у) и границы Е„(х,у), для которых величины to , успр, рассчитанные численными методами, совпадают с экспериментальными значениями ^, усвр .
Предположим, что мы умеем приближенно решать обратную задачу, то есть умеем строить (численно или аналитически) такой оператор (р, что £ ~ <р(йО, откуда (р -1, где /-1 - оператор, обратный оператору /. Ранее такой оператор был построен. Он является приближенным, так как не учитывает кривизну луча в оцениваемом слое и отличие интегральной скорости успр от дифференциальной уе, которая соответствует интегральной при нулевой базе (Пузырев, 1989). Другими словами, оператор построен в предположении, что значения интегральной скорости успр определены при малых базах годографов ОГТ, а также в предположении локальной однородности оцениваемого слоя.
Рассмотрим несколько подробнее понятие интегральной скорости усор и связанное с ним понятие дифференциальной скорости. Предположим, что шаг А1 = 1/+1 -1/ достаточно мал, так что сумму по / в (14) можно заменить на интеграл (на практике это условие обычно выполняется) (Мешбей, 1985). В этом случае скорость уср,р, минимизирующая функцию Ф, зависит от Ь1, Ь2, где Ь1, Ь2 - минимальное и максимальное значение расстояний I/ взрыв-прием. Это значит, что усвр = усвр(^1, Ь2). Если Ь1, Ь2 стремятся к нулю, то усор(^1, Ь2) стремится к некоторому значению, которое будем обозначать через Уе, и которое называется предельной эффективной скоростью. Если годограф ОГТ является гиперболическим (что соответствует однородной среде с плоской отражающей границей), то Ущр не зависит от Ь1, Ь2, и, следовательно, усвр = Уе. В более сложных моделях реальных сред (в том числе и горизонтально-слоистой модели с однородными слоями) усор отличается от уе.
Уменьшение базы годографа при проведении скоростного анализа приводит к сильному уменьшению устойчивости определения скорости, но, с другой стороны, при этом в большинстве случаев годограф ОГТ близок к гиперболе, поэтому отличие интегральной скорости от дифференциальной не очень велико. Влияние неучета кривизны лучей в самом оцениваемом слое при определении его скорости и мощности также невелико (Бляс, 1991). Это значит, что неучет отличия интегральной скорости успр от дифференциальной уе и неучет кривизны луча в оцениваемом (!) слое не приводит к большим погрешностям в определении пластовой скорости и мощности, поэтому оператор в самом деле дает приближенное решение обратной задачи, то есть <р~/-1. С другой стороны, если для решения обратной задачи применять только оператор (р, то получим систематические погрешности, которые хотя и невелики, но могут существенно повлиять на результаты геологической интерпретации, особенно в районах с платформенным залеганием пород и малоамплитудными структурами.
Для учета отличия усор от Уе рассмотрим следующий итерационный алгоритм:
£т+1 = £т - <р(/(Щ) + ), £1 = т = 1,2,... . (15)
Действие данного алгоритма сводится к следующему. На первом шаге оператор ф применяется к экспериментальным данным ё (к значениям интегральной скорости успр и времени to по профилям ОГТ), его применение дает приближенное решение обратной задачи - получаем Б1= (р(й). Так как оператор ф дает точное решение только для предельной эффективной скорости и локально-однородного оцениваемого слоя, то полученное значение £1 будет приближенным. К полученной на первом шаге модели £1 применяется оператор / решения прямой задачи - для данной модели численно рассчитываются годографы ОГТ (спр( I), которые аппроксимируются гиперболами. Параметры гиперболы рассчитываются из условия минимума функции (14). К полученным таким образом значениям (численно рассчитанным для модели £1, /02, усор) применяется оператор (р, который дает модель £1 = <Р(/(£1)). Затем модель корректируется добавлением величины АБ1 = £/: Б2 = Б1 + А81. На том шаге находим 8т' = (р(/(£т)), А£т = Бт- 8т' и новую модель Бт+1 = Бт+ А£т.
Если (р = /(то есть на первом шаге получаем точное решение £1 = Я* обратной задачи), то £1' = <Р(/(£1)) = £* и ^£1 = £1- £1 = 0. В этом случае все модели Ят совпадают, и первая итерация дает точное решение обратной задачи. Если <р^/-1, то на каждом шаге алгоритма получаем различные модели
£2,..., «^иь... .
Возникает вопрос о сходимости данного итерационного алгоритма и о модели, к которой он сходится. Исследуем сходимость алгоритма. Равенство (15) можно записать в виде
Бт+1 = А(Бт), А(Б) = Б - <р(/(Б)) + <р(й), т = 1,2,.. . (16)
Найдем производную оператора А. Из равенства (16) следует, что А' = Е - р' /', где Е -единичный оператор. Так как (р то р' ')'1, откуда Е ~ (р'/' и \\А'\\ = ||Е - р'/'|| « 0, где - норма оператора А. Отсюда следует, что оператор А является сжимающим (его норма меньше единицы), и, по теореме об операторе сжатия, итерационный алгоритм (16) сходится (Колмогоров, Фомин, 1968). Скорость сходимости зависит от ЦА'Ц и определяется формулой
||Бт - Б*|| < РГ ||А(Б1) - й|| (1- |И'|| ) -1 ,
где Б* - модель среды, к которой сходится последовательность Бп при п Переходя в равенстве (15) к пределу при п ^ да, получим, что Б* является решением уравнения
(р( / (Б)) = <р(й).
Если оператор р имеет обратный, то из последнего равенства следует, что Б* удовлетворяет уравнению (13), то есть Б* является решением поставленной задачи и, следовательно, не зависит от выбора оператора (р. Отсюда следует важный вывод: оператор р можно конструировать в более простой модели среды (или более простой схеме наблюдений) в предположении, что нам известны исходные данные, которые на самом деле по реальным сейсмограммам не могут быть определены. Для получения решения обратной задачи (решения уравнения (13)) необходимо в операторе /решения прямой задачи учитывать более сложную модель, чем та, которая рассматривается при построении оператора (р. В рассматриваемой задаче это означает, что при решении обратной задачи (построении оператора (р) можно использовать модель с локально-однородным оцениваемым слоем (для которого только и можно построить конструктивный алгоритм определения скорости и мощности) и считать, что известна предельная эффективная скорость (хотя на практике она не может быть устойчиво определена по годографу, так как ее несмещенная оценка имеет бесконечную дисперсию). При расчете же траекторий и времен в прямой задаче необходимо учитывать кривизну лучей, а затем рассчитывать интегральную скорость по гиперболической аппроксимации годографов ОГТ.
Итак, данный алгоритм позволяет определить скоростную модель с учетом кривизны лучей в горизонтально-неоднородных слоях несмотря на то, что инверсный метод предполагает локальную однородность оцениваемого слоя.
5. Оптимизационный алгоритм решения обратной задачи
Рассмотрим применение методов оптимизации к определению скоростной модели среды по годографам отраженных волн. Так как на практике при поверхностных наблюдениях в большинстве случаев отраженные волны выделяются с помощью скоростного анализа, а не их потрассной корреляцией, то нам известны не сами годографы, а их аппроксимации, чаще всего гиперболические (Мешбей, 1985). В этом случае по выделенным гиперболам можно рассчитать годографы и считать, что нам известны времена прихода отраженных волн, в то время как при скважинных сейсмических наблюдениях обычно известны сами годографы.
Пусть, как и ранее, Рк(х,у), Ук(х,у,1) - функции, описывающие границы слоев и пластовые скорости среды, скоростную модель которой обозначим через Б. Через / обозначим номер точки наблюдения, /=1,2,..., I; р - номер волны в множестве выделенных волн, р=1,2,...,Р. Под /-ой точкой наблюдения понимается пара источник-приемник, имеющая номер / в общей нумерации. Через ир обозначим известное из экспериментальных данных время пробега р-ой волны в /-ой точке наблюдения, иР(Б') - "теоретическое" время пробега волны, рассчитанное для модели Б'. Будем рассматривать трехмерные модели сред, в которых границы слоев Рк(х,у) и пластовые медленности У^_1(х,у,г) описываются линейными комбинациями заранее выбранных базисных функций:
Ук~\х,у,2)= «у(к) р(х,у,2),
к=1,...,К (17)
Рк(х,у)= щ(х,у,х).
Метод оптимизации сводит обратную задачу к нахождению коэффициентов а^, р/к) из условия минимума функционала Ф, определенного равенством
Ф(а,Р) = I (tp- tip(a,fi))2qip, (18)
где qip - веса, характеризующие погрешность определения tip, а, Р - векторы с координатами, состоящими из коэффициентов разложений (17):
сс = (о*1 1, СС2 1,... 1, сс/\..., ccjJ^,..., а/к>, ccjccJJ
(19)
Р=(Р/1}, рР,.., j(1), рЯ., j(2),..., рЯ p2(K>,.., J<K>).
Через Pn(£n(,'p>, riJ'-p>, F„(^n(',p, r]„(',p)) обозначим точки пересечения траектории p-ой волны с границами слоев, перенумерованных в порядке встречи с лучом, n=1,...,N. Время tip(a,P) пробега p-ой волны в /-ой точке наблюдения в модели можно представить в виде
N
T,p(a,Р) = I tn($Jv>, ^Я, &hp), чЯ) = T(ip)(a,p, ?p), (20)
n=1
где tn(^n-1(l,p>, rnJ'-p>, ^Я, t7n(lp>) - время пробега волны между точками Pn-1, Pn, то есть время пробега волны в n-ом слое, ¿*,р> = (4я, ^я,..., 4n'p>), ^hp>=(^1(bp>, ri2(hp),..., ^я) - векторы, описывающие траекторию p-ой волны в /-ой точке наблюдения. В соответствии с принципом стационарности
m<hp>/d4n(hp>=0, dt^/dn™ = 0, n=1,...,N. (21)
Итерационный алгоритм нахождения минимума функции Ф, которая описывает квадратичное отклонение теоретических времен от реальных, сходится за небольшое количество итераций. На каждой итерации выполняется линеаризация функционала Ф и рассчитываются траектории (число расчетов траекторий совпадает с числом итераций), для нахождения минимума линеаризованного функционала можно применить методы квадратичной минимизации, которые сходятся за конечное число шагов.
Рассмотрим да-ый шаг алгоритма минимизации функционала Ф. Пусть S0 = (а0,р0) - модель, полученная на да-ой итерации. Через {1;Я, г/Я} обозначим траекторию p-ой волны в модели S. Разложим функцию Tip(a, Р, £,р> , rf,p>) в ряд Тэйлора в точке (а0,р0, £о(,р>,цЯ). С учетом того, что ¿¡=£(а,Р), ri=ri(a,P) и равенств (21) получим, что
T(,p) = Toip> + жЯ/да(а-ао) + ^Я/др(Р-Ро) + R(a, (22)
где R(a,P) - величина второго порядка малости по отношению к ||а-ао||, \\Р-Ро\\; dT0( lp/da, ¿T0( lpp/dp -векторы с координатами
JT0(,р)/да = (JT0(,'p)/da1(1), JT0(^/даЯ..., JT0(,р)/даЯ, dT0(/р)/да1(2)...., JT0(,p)/daju(1}), JT0(,p)/dp = {0TO(hp>/dj31(1>, dT0(^/дрЯ,.., dT0(,р)/дрЯ, dT0(,p)/dp1(2)..., dT0(,p)/Щи(1)),
"0" в T(,'p и ее производными означает, что они берутся в точке (а0,Р&1;Я, 'hp).
Подставим (22) в (18) и ограничимся учетом линейных членов по а-а0, Р-Ро. Тогда получим, что
Ф=Т [tp- T<ip>- dT0(гр)/да(а-00) - T(lp)/dp(fi-р0)\\гГ= Ф(*Р). (23)
'p
Так как функция T,р> и ее производные сГ0(,р>/.да, сГ0(lpp/dp вычисляются в известной модели S0, то неизвестные а, Р входят в Ф квадратично. Это значит, что задача минимизации функции Ф сводится к минимизации квадратичной функции по координатам векторов а-а0, Р-Р0. Эта задача решается за конечное число шагов, например, методом сопряженных градиентов.
Пусть Atip = tip-T0p - разность между наблюдаемыми временами tip и временами рассчитанными для модели S0, Аа =а - а0, АР=Р - р0 - искомые добавки к значениям а, Р модели S0, полученной на предыдущей итерации. Тогда функционал Ф можно записать в виде
Ф = Z(At,p- JT0(l'p)/daAa-Ж0( 1р)/дрAP)2qhF hp
В Ф входят производные функции по коэффициентам aj^, f¡jkk разложений пластовых скоростей и границ слоев. Для '-ой точки наблюдения эти производные считаются вдоль траектории, соединяющей источник и приемник, k - номер слоя. Для нахождения этих производных достаточно решить один раз прямую задачу - найти траектории в модели S0 = (a0,fi0), а затем явным дифференцированием функции (20) найти производные времени Т(,"рР по коэффициентам afkk, fijkk разложений (17) границ и скоростей.
В функцию Т(,"р входят времена tk пробега волны во всех пересекаемых слоях, k = 1,2,...,n. При дифференцировании по коэффициентам ajkk скорости k-ro слоя останутся только времена пробега в этом слое, если же мы находим производную по коэффициентам J3jkk n-ой границы, то остаются производные времен в слоях, разделяемых этой границей, то есть в слоях с номерами k-1 и k.
Если слои можно считать локально-однородными, то время tk в k-ом слое можно найти по явной формуле
tk = {[Fk(4 Щ) - Fk-1(&1, %-1)]2 + (&-&-1)2 + (%-%-1)2}2mk, (24)
где mk=Vk-l((^k-1 + ^k)/2,(^k-1 + 4k)/2,(Fk-1+Fk)/2) - значение медленности в средней точке отрезка [Ры, Р]. Подставляя в (24) равенства (17), получим явную зависимость времени tk от коэффициентов otfk>, Pjk~11>, Pjkk по которой легко найти производные. Для среды со слабо неоднородными слоями в Приложении 2 получена явная формула для времени пробега волны в слое, по которой также можно найти производные, если в нее подставить равенства (17). Если же слой существенно неоднородный, и изменение скорости в нем нельзя учесть методом возмущений, то для нахождения траекторий надо применять численные методы решения обыкновенных дифференциальных уравнений, а для нахождения производных можно применить метод, предложенный автором (Бляс, 1985).
Таким образом, нахождение минимума функционала Ф сводится к нахождению минимума квадратичной функции Ф на каждом шаге алгоритма. Для нахождения коэффициентов этой квадратичной функции необходимо один раз решить прямую задачу - рассчитать траектории волн для заданной схемы наблюдений, после чего найти производные функции T0( 1,рР, входящие в (23). Нахождение этих производных сводится к нахождению производных времен пробега в отдельных слоях. Численные расчеты показывают, что если пластовые скорости модели нулевого приближения (априорной модели) отличаются от истинных скоростей не более чем на 15 - 20%, то алгоритм сходится за 3 - 5 итераций.
Приложение 1. Эффективные скорости, определяемые по кинематическим и динамическим годографам
В теоретических исследованиях, касающихся интегральной скорости vCDP, обычно предполагается, что эта скорость найдена по гиперболической аппроксимации годографа tCDP(l) ОГТ по методу наименьших квадратов (Пузырев, 1979; Гольдин, 1993). Другими словами, vCDP определяется как значение параметра v, при котором функция Ф, определенная равенством
N
ф&) = I [ </,) - (0 + /'W2]2, (П1.1)
,=1
достигает своего наименьшего значения (l, - удаление взрыв-прием на '-ом приемнике). В практике сейсморазведочных работ методом отраженных волн сами годографы ОГТ обычно неизвестны, и vCDP определяется в результате скоростного анализа сейсмограмм ОГТ (Мешбей, 1985; Урупов, Левин, 1985). Метод скоростного анализа состоит в накапливании (в том или ином виде) в некотором скользящем
временном окне сейсмических трасс u,(t) с некоторыми перебираемыми подвижками в, (i - номер трассы
2 2 2 1/2
сейсмограммы ОГТ). Обычно 0j = (t0 + /vCDP ) - t0, но иногда, в случае сложного строения покрывающей толщи, когда необходимо учесть отличие годографа ОГТ от гиперболы, применяются более сложные аппроксимации, зависящие от трех перебираемых параметров (Гольдин, 1993). Результат накапливания представляется в числовом виде, то есть для фиксированного временного окна каждому набору задержек (в случае гиперболической аппроксимации - каждому значению, определяющему эти задержки) ставится в соответствие некоторое число E (рис.2б), дающее оценку когерентности (синфазности) колебаний вдоль кривой (рис.2а). Оператор, переводящий результат накапливания (или некоторой процедуры оценки когерентности трасс в заданном окне, включающей в той или иной форме суммирование трасс) в число, называется оператором регулируемого направленного анализа (РНА). Значения оператора РНА, полученные для некоторого множества значений vCDP, называются скоростным спектром, рис.2б (Урупов, Левин, 1985; Маловичко, 1990; Карасик, 1993). Так как оператор скоростного
анализа зависит также от окна суммирования, то величина Е зависит от двух переменных: от скорости суммирования и времени /0 - центра окна суммирования (рис.2в).
Если в окне анализа сейсмограммы имеется одна регулярная волна с гиперболическим или
2 2 2 1/2
близким к гиперболическому годографом, то при совпадении кривой суммирования (/0 + /усор ) с годографом этой волны происходит синфазное накапливание, и оператор РИА принимает максимальное значение. Этому максимальному значению соответствует некоторое значение успр, получаемое в результате скоростного анализа. Это значение усвр в дальнейшем используется как для получения суммарных разрезов ОГТ, так и для определения пластовых скоростей в результате решения обратной кинематической задачи тем или иным методом.
а
///у/
-> У ОГТ
Уогт
Рис.2. Скоростной анализ сейсмограммы.
Несмотря на большое количество различных операторов скоростного анализа, оценки уСпр, получаемые этими операторами, близки друг к другу, а при отсутствии помех (как регулярных, так и случайных шумов) вообще совпадают. Фактически все эти операторы являются операторами второго порядка и выражаются через одни и те же величины, связанные с суммой трасс в заданном окне. Подробно эти вопросы рассмотрены во многих работах, обзор этих работ и их анализ дан Маловичко (Маловичко, 1993). Здесь мы рассмотрим энергетический оператор Е, дающий энергию суммотрассы. Аналогично можно рассмотреть действие других операторов.
2 2 2 1/2
Пусть и() - /-ая трасса сейсмограммы ОГТ, 1(у)=((0 + 1/ /уСор ) - кривая суммирования. Вместо коэффициента V удобнее рассматривать параметр Ь = 1/у2, входящий в подкоренное выражение линейно, в то время как сама скорость суммирования входит в подкоренное выражение нелинейно. Кроме того, введение параметра Ь позволяет рассматривать и так называемый случай переворачивания годографов ОГТ, вызванного резкими криволинейными преломляющими границами в покрывающей среде (Бляс, 1991). Для параметра Ь кривая суммирования т(Ь) имеет вид т(Ь) = (/02 + Ы,2)1/2, и энергетический оператор Е скоростного анализа дается формулой
Т N
Е(Ь) = I [ I и (/ + (0 + Ы'ГШ (П1.2)
0 /=1
Здесь для удобства исследований предполагается, что шаг временной дискретизации Л/ достаточно мал, и вместо суммы по / при нахождении энергии суммотрассы берется интеграл. На практике это условие в подавляющем большинстве случаев выполняется, так как работы методом ОГТ проводятся приемными устройствами с шагом дискретизации 1 мс, что составляет сотые доли от периода полезной волны. При интегрировании будем считать, что полезная волна полностью попадает в интервал интегрирования.
Предположим, что в окно анализа попадают две регулярные интерферирующие волны, то есть,
u,(t) = ft - Tf (/,)) + g(t - тг(/,)), /=1,2,...,«. (П1.3)
Заметим, что в принципе можно рассмотреть произвольное число интерферирующих волн, а также учесть наличие нерегулярного (некоррелируемого по трассам) шума, если под слагаемым g понимать не одну волну, а сумму нескольких регулярных волн и случайного шума. Для сокращения записи мы ограничиваемся учетом интерференции двух регулярных волн.
Пусть f (t) - сигнал полезной однократно отраженной волны, чей годограф описывается функцией
T(l) = [со + C112 + C214] 1/2
(П1.4)
Здесь предполагается, что реальный годограф описывается равенством (П1.4), то есть более сложной кривой, чем гипербола. Автором (Бляс, 1988) показано, что во многих случаях годограф ОГТ с высокой степенью точности можно описать формулой (П1.4). Так как мы исследуем результаты гиперболического скоростного анализа, то здесь уместно заметить, что если годограф ОГТ сильно отличается от гиперболы, то в операторе РНА необходимо учитывать это отличие, иначе максимумы спектров скоростей будут существенно искажены, и, следовательно, скорости суммирования невозможно будет использовать в дальнейшем как для получения скоростной модели среды, так и для суммирования ОГТ.
Подставив (П1.3) в (П1.2), получим формулу, дающую явную зависимость оператора Е от параметра Ь:
Т N
Е(Ь) = [/0 - т/(I) + (0 + ы2)1/2} + Ж - гД) + (0 + Ь/2)1/2)]2Н (П1.5)
0 1=1
Нас интересует точка максимума этой функции, то есть такое значение Ь* параметра Ь, в которой ёЕ/ёЬ = 0. Подставляя в это равенство представление (П1.5), получаем следующее уравнение, которому удовлетворяет точка Ь* максимума оператораЕ:
Т N N
1{[ I& - гД) + (0 + ы2)1/2)+1 Ж - гД) + (0 + Ь1?Г) ] X
0 = = (П1.6)
N N
х [ If (- т(1) + (0 + W2)1/2)+ I g '(t- Tg(ld + (0 + bl2)1/2} ] } [l2/(to2 + bl2)
i=1 i=1
Точное аналитическое решение данного уравнения получить не удается, но, используя метод возмущений, можно получить явное решение с любой точностью. Анализ данного решения позволяет исследовать отличие уСвр, найденного по годографу ОГТ методом наименьших квадратов (то есть чисто кинематическое значение интегральной эффективной скорости, полученное при минимизации функции Ф), от значения уСпр, полученного при скоростном анализе сейсмограмм. Интегральную скорость, определенную путем минимизации функции Ф, будем называть кинематической эффективной скоростью и обозначать ее через ^Свр, формула для кинематической эффективной скорости будет получена ниже. Эффективную скорость, найденную при скоростном анализе сейсмограммы ОГТ, будем называть динамической скоростью (или скоростью суммирования) и по-прежнему обозначать уСвр. Очевидно, что кинематическая скорость wCDP зависит только от годографа ОГТ отраженной волны (если он отличен от гиперболы) и от базы наблюдения (от промежутка аппроксимации реального годографа гиперболой), в то время как динамическая скорость у^р зависит также от динамических характеристик волн - формы их сигналов/(/), g(t) и от характера интерференции.
Применение метода возмущений накладывает некоторые ограничения на модель сейсмограммы, а именно - рассматриваемая модель должна "не очень сильно" отличаться от базовой модели (модели нулевого приближения), для которой решение рассматриваемого уравнения легко находится. Это приводит к следующим предположениям относительно модели сейсмограммы, которые на практике часто выполняются, по крайней мере, для хорошо выраженных отражений и платформенных условий залегания пород. Во-первых, предполагается, что выполняется с214 « 10+с112, то есть влияние третьего коэффициента с2 на годограф ОГТ намного меньше, чем влияние суммы первых двух. Это условие на практике в большинстве случаев выполняется, по крайней мере при удалениях, не превышающих глубины отражающей границы (Бляс, 1988). Строго говоря, понятие малости отклонения годографа от его гиперболической аппроксимации зависит от преобладающей частоты сигнала, само отклонение
должно выражаться в долях периода сигнала, подробнее этот вопрос будет рассмотрен в конце приложения. Во-вторых, предполагается, что при Ь=с1, то есть при суммировании по кривой, близкой к годографу полезной волны, энергия суммарного сигнала больше энергии помехи. Другими словами, предполагается, что суммирование по гиперболе, аппроксимирующей годограф полезной волны, позволяет выделить эту волну на фоне помехи, то есть разрастание на скоростном спектре вызвано полезной волной. Численные расчеты на моделях показывают, что, если отличие годографов Т/ и тг на последних каналах (при максимальном расстоянии взрыв-прием) достигает величины порядка периода сигнала, то при равных амплитудах это условие выполняется.
Согласно методу возмущений, малость отклонений годографа ОГТ Т/ от гиперболы и малость суммарной волны помехи (малость помехи по сравнению с сигналом на суммотрассе) перенесем в малый безразмерный параметр е, для чего вместе с исходным уравнением (П1.6) рассмотрим вспомогательное уравнение
Т N N
1{[!/(г - т(/„£) + (0 + Ь/2)1/2} + £Т & - гД) + (0 + Ь/2)1/2} ] х 0 '=' = (П1.7)
N N
х [I/'0- т/(1,ё) + (0 + Ы?Г) + еЪ- т,(/) + (0 + Ь/2)1/2}]} [/2/(/02 + Ь/2)] Лг,
1=1 1=1
где т/у/ьё)=[оо + С1/2 + £С2/4]1/2.
Прежде всего покажем, что при е=0, то есть при отсутствии волны-помехи и гиперболическом годографе полезной волны, суммирование по годографу волны приводит к максимуму энергетического оператора, то есть к максимуму спектра скоростей. Во временной области этот факт обосновывается с помощью неравенства Коши - Буняковского (Колмогоров, Фомин, 1968): для любых чисел хг, уг, г=1,2,...,п, выполняется неравенство
п п п
( Е Xi Уг )2 <Е xi Е Уi ,
1=1 г=1 г=1
причем равенство достигается только при условии, что хг= уг для всех г=1,2,...,п. Подставляя в данное равенство xi = /($ - $), и учитывая, что I/2(/ - $) Ж не зависит от г, получим, что
Т N Т
I (I/( г - в'))2 Л < N I/2(0^,
о г=1 о
где в'== (г°+Ь/2)112- т/(/г,0). Строгое равенство выполняется в случае, когда/(г-в) принимают одинаковые значения для всех г, что соответствует синфазному суммированию.
Итак, при е =0 уравнение (П1.7) имеет решение Ь0=с1, поэтому согласно методу возмущений можно найти решение этого уравнения в виде ряда по степеням £:
Ь0 = с1 + еЬ1 + £2Ь2 + ... . (П1.8)
Для нахождения коэффициентов Ьг, г=1,2,... ряда (П1.8) подставим его в уравнение (П1.7) и разложим левую часть в ряд по степеням £ Так как правая часть получившегося тождества не зависит от £, то все коэффициенты при степенях е в левой части также должны равняться нулю. Приравнивая их к нулю и учитывая, что с1 - решение уравнения (П1.7) при £=0, получим линейные уравнения для нахождения коэффициентов Ьг, г=2,3,... . В каждое такое уравнение, получающееся приравниванием к нулю коэффициента при ё, входят коэффициенты Ь1, Ь2, ..., Ьг-1, и линейно Ьг. Отсюда следует, что решая последовательно получающиеся уравнения, можно найти сколько угодно коэффициентов ряда (П1.8). Если считать, что сигналы описываются аналитическими функциями, то по теореме об аналитичности неявной функции (ТФКП) ряд (П1.8) также описывает аналитическую функцию в некоторой окрестности нуля, то есть является сходящимся рядом. Хотя в принципе можно найти сколько угодно членов ряда (П1.8), для аналитических исследований достаточно ограничиться первыми двумя слагаемыми. Численные расчеты показывают, что они дают значение эффективной скорости с достаточно высокой точностью.
Ограничиваясь учетом линейных членов по £И переходя от вспомогательного уравнения (П1.7) к исходному уравнению (П1.6) (которое не зависит от е), получим, что энергетический оператор Е достигает максимума при значении ЬСвр, которое дается равенством
N N N N N
Ьспр = С,1 + [( I 12 /D}l2 )2 - N Т 14 /D1 ]-1 { [ !(Л4 / DI1/2) (112 /DI1/2) - N (I Ц /D1)] С2 1=1 1=1 1=1 1=1 1=1
Т Т N Т
- [2 / \ /'2(0 Ж] [ \ /'(О О(Г) Ж I I,2 /D}l2 + N У) ем Ж ]},
0 0 ,=1 о
где
Di = со + С121 \ О(Г) = I gi (Г - + (со + С1 12)1/2},
01 (О = I g'i (( - + (со2 + с^Г) I2 /D}'2.
Для упрощения данного равенства воспользуемся следующим приемом, позволяющим перейти от сумм к интегралам. Предположим, что шаг Л1 = 1+ - I, между трассами достаточно мал, и суммы по , можно заменить на интегралы в пределах базы скоростного анализа, например,
N Ь
I [/2/(со + с /2)] * (М/Ь) I [/2ё// (Со + с/)] .
,=1 о
На практике это условие в подавляющем большинстве случаев выполняется, так как в настоящее время работы методом ОГТ проводятся с приемными устройствами, в которых расстояние между приемниками составляет сотые и даже тысячные доли от длины расстановки.
Далее, введем в рассмотрение безразмерную величину уЬ по формуле уЬ = с1Ь2/ со. Для горизонтально-слоистой среды с мощностями и пластовыми скоростями /Ь=(1/4) (Ь /Иг)2, где
N N
Ие = ( I ИкУк I кк /Ук )1/2 к=1 к=1
- эффективная глубина, близкая к глубине Н отражающей границы. Для сред с небольшими углами наклона (соответствующими платформенным условиям осадконакопления) величина уЬ также примерно равна одной четвертой квадрата отношения длины расстановки к глубине отражающей границы, формирующей полезную волну, то есть много меньше единицы.
Беря интегралы в явном виде, разлагая получающиеся выражения в ряд по степеням у и ограничиваясь учетом первой степени у, равенство (П1.9) можно записать в более простом виде:
Т
ЬС№ = с1 + (6 /7)с2Ь2 + (45 с0 /2Ж4) (1+ бу/7) / [ \/'2(0 Ж] х
0 (П1.11)
Т Т
X [Ь2 (1 - 3уь/10) / (30 \/'(0 0(0 Ж + У) от ж ].
0 0
Данная формула дает явное (приближенное) представление точки максимума скоростного спектра в случае небольших отличий годографа ОГТ от гиперболы и небольшой энергии волн-помех на суммотрассе ОГТ. Первое слагаемое дает значение параметра Ь, совпадающее с коэффициентом с1, то есть дифференциальным значением, получающимся на бесконечно малой базе анализа при гиперболическом годографе и отсутствии помех. Второе и третье слагаемые в правой части формулы (П1.11) дают смещение оценки коэффициента Ьспр относительно с1=1/уг2. Второе слагаемое описывает смещение, вызванное отличием годографа полезной волны от гиперболы, третье дает влияние волны-помехи, которая искажает оценку эффективной скорости.
Если волна-помеха отсутствует, то из (П1.11) получаем
ЬcDP = 1 / у^р = С1 + (6 / 7 )Ь2С2. (П1.12)
Прежде всего отметим, что при отсутствии помехи величина Ьспр, а вместе с ней и УтР=(1/с_;)1/2 не зависит от формы сигнала, а только от базы скоростного анализа и коэффициента с2, описывающего отличие годографа полезной волны от гиперболы.
Для анализа отличия кинематической и динамической эффективных скоростей найдем значение коэффициента Ь, минимизирующего функцию Ф, определенную равенством (П1.1). Как и в случае динамической эффективной скорости, будем считать, что годограф т(1) описывается равенством (П1.4).
(П1.9)
(П1.10)
Для получения явного представления коэффициента Ь также применим метод возмущений, для чего вместе с исходной функцией Ф рассмотрим вспомогательную функцию Ф, зависящую от малого безразмерного параметра £ и определенную равенством
N
Ф(У) = Е [(Со + Ы2 + £€214 )1/2 - (Со + (П1.13)
1=1
Точка минимума данной функции удовлетворяет уравнению
N N
сФ/дЬ = Е [(с0 + с112 + с214)1/2 /(с0 + С/2)1/2]/2 = Е I2. (П1.14)
г=1 г=1
При £= 0 оно имеет очевидное решение С = с1, поэтому будем искать решение данного уравнения в виде ряда по степеням £:
С = С1 + £ С1 + £2С2 + ... .
Подставим данный ряд в уравнение (П1.14), разложим левую часть в ряд по степеням £ и приравняем коэффициенты при степенях £ к нулю. Получим уравнения, решая которые, можно последовательно найти сколько угодно коэффициентов ряда (П1.8). Ограничиваясь учетом первой степени £, получим
N N
Ссор = С1 + С2 Е ^/(со + С112) / Е ///(со + С112).
1=1 1=1
Как и ранее, для упрощения данного равенства заменим суммы на интегралы и ограничимся учетом первой степени параметра у. Тогда данная формула упрощается и принимает вид
ССОР = 1 / ЯВР = С1 + (5 / 7)^. (П1.15)
Из равенств (П1.15) и (П1.12) следует, что кинематическая и динамическая скорости даже в отсутствие помех отличаются, но это отличие невелико, и легко может быть учтено при теоретических исследованиях - для этого достаточно во всех формулах, куда входит интегральное значение скорости, вместо Ь подставить (6 / 5)1/2 Ь .
Теперь рассмотрим важный для практики вопрос о "малости" коэффициента с2, при которой можно воспользоваться разложением в ряд по степеням £ в уравнении (П1.7) или, что то же самое, в ряд по степеням с2. Качественно ясно, что понятие малости с2, или, что то же самое, малости отличия годографа полезной волны от гиперболы, зависит от частоты сигнала - чем выше его частота, тем сильнее влияние негиперболичности годографа на суммарный сигнал. Для получения количественных оценок рассмотрим функцию
КО = Е/(г - (Со + С11 2 + С214)1/2 + (Со + С112)1/2)
в спектральной области. Пусть - преобразование Фурье функции/(1), то есть
да
Р(®) = (/(0 егМ Ж.
о
Тогда из свойств преобразования Фурье следует, что
= ад Е е-Шч,
2 1/2 2 4 1/2 2
где - спектр функции ^(0, тк = (со+ с11) - (со+ с11 + с2Г) . С точностью до 0(с2 ) имеем
равенство
тк = - ( 1 / 2 )С214 / (со + С112)1/2,
поэтому разложение в ряд по степеням с2 с точностью до с2 равносильно разложению е в ряд по степеням тс точностью до (га>тк)2. Отсюда следует, что возможность учета только первой степени с2 (с отбрасыванием более высоких степеней) фактически соответствует возможности отбрасывания степеней
(/®т) выше первой при разложении функции е,т в ряд по степеням 1ог. Для функции ет справедливо разложение
ег = 1 + 2/1! + 22/2\ + 23/3\ + ... + 2п/п\ + ... .
Можно считать, что е2« 1+2 при | х | < 0,5, то есть в разложении функции еШТ в ряд по степеням га>т можно ограничиться учетом только первой степени 1а>т, если выполняется неравенство | сот \ < 0,5.
Далее, пусть а>„ - преобладающая частота сигнала, Тп - период, соответствующий этой частоте, то есть Тп = 2п/а>. Тогда условие возможности учета только первой степени с2 при разложении в ряд принимает вид 2пт/Т„< 0,5, или приближенно т< Г/10. Это значит, что применение метода возмущений для описания отличия годографа т(Г) от гиперболы при анализе спектра скоростей возможно, если это отличие не превосходит 1/10 видимого периода. В практике морских сейсморазведочных работ преобладающий период полезных (однократно отраженных продольных) волн для глубин порядка 2 - 3 км составляет 30 - 40 мс, поэтому отличие годографа от гиперболы не должно превосходить 3 - 4 мс, что, как показывают численные расчеты, выполняется для большинства нефтеносных районов (Пузырев, 1979; Маловичко, 1990).
Приложение 2. Уравнение траектории и время пробега волны в неоднородном слое
Пусть в неоднородном слое со скоростью у=у(х,у,х) заданы две точки Q1(41,r|1,^1) и Q2(42,r|2,^2). Требуется найти траекторию у. х=х(х), у=у(х), соединяющую эти точки, и время пробега волны. Сначала получим приближенное представление траектории, а затем найдем время в виде ряда по степеням малого безразмерного параметра. Функции х=х(х), у=у(х), описывающие траекторию у, удовлетворяют уравнениям Эйлера
дЕ/Зс - ё(сЕ/Зс')/ёх = о, дЕ/су - ^сЕ/су')/ёх = о, (П2.1)
где
Е(х,у,х,х',у0 = (1 + х'2 + у'2)1/2 / у(х,ух) , (П2.2)
и функции х(х), у(х), удовлетворяют краевым условиям
х(£) = 41, у(£) = ГЦ. (П2.3)
Подставляя (П2.2) в (П2.1) и раскрывая производные с учетом того, что х и у являются функциями 2, получим, что функции х(х), у(х) удовлетворяют системе дифференциальных уравнений второго порядка
(^2Е /Зс'2)х" + (д2Е/сх'суг)у"+ (д2Е/сх'сх) х'+ (д2Е/сх'ду)у'+ (д2Е/&'&) - дЕ/сх = о ,
(П2.4)
(д2Е/дх'ду')х" + (д2Е/ду'2)у" + (д2Е/ду'дх) х'+ (д2Е/ду'ду)у'+ (д2Е/су'дх) - дЕ/ду = о . Для однородной среды, когда скорость V не зависит от х, у, 2, система (П2.4) переходит в уравнения
(д2Е/дх'2)х" + (д2Е/ск'ду)у" = о ,
(д2Е/сх'су')х"+ (д2Е/ду'2)у" = о ,
откуда х"(г)=0, у"(г)=0. Тогда с учетом краевых условий (П2.3) получаем, что
(П2.5)
у(2)=ГЦ + [(^1)/(£-£)](2-а
и траектория ^представляет собой прямую. Если слой неоднородный, то точное решение системы (П2.4) удается найти только для среды с линейным изменением скорости (Пузырев, 1979). В то же время, как показано в (Бляс, 1991), в верхней части разреза основную роль играют нелинейные (квадратичные) изменения пластовых скоростей, поэтому важно иметь явное представление времени пробега волны в слоях с нелинейным изменением скоростей.
Пусть h(x,y,z) - медленность, то есть величина, обратная скорости: h(x,y,z) = 1 / v(x,y,z). Тогда уравнения (П2.4) можно переписать в виде
nx''=(1+x' 2+y' 2){ch/cX - ch/dz),
(П2.6)
ny"=(1+x' 2+y' 2)(ch/cy - ch/cZ).
Для получения явной формулы для траектории волны и ее времени пробега вместе с исходной моделью медленности, описываемой функцией n(x,y,z), рассмотрим вспомогательную модель с медленностью, зависящей от малого безразмерного параметра S:
h(x,y,z,s) = h0 + sm(x,y,z). (П2.7)
Для модели с медленностью h(x,y,z,s) система (П2.7) принимает вид
(n0 + sm(x,y,z)) x" = s(1+x' 2+y' 2)(3n/Sx - 3n/dz),
(П2.8)
(h0 + sm(x,y,z)) y" = s(1+x' 2+y' 2)(3n/3y - 3n/dz).
Согласно методу возмущений, решение краевой задачи (П2.3), (П2.8) ищем в виде ряда по степеням г:
x (z,s) = x0(z) + sx1(z) + s2x2(z) + ... ,
(П2.9)
y(z,s) =yo(z) + syi(z) + S2y2(z) + ... .
Для нахождения коэффициентов x(z), y(z) рядов (П2.9) подставим их в уравнения (П2.8), разложим левые и правые части в ряды по степеням s и приравняем коэффициенты в левой и правой частях при одинаковых степенях е. Получим дифференциальные уравнения относительно x(z), y(z), ¿=0,1,2,... . Для х0, у0 получим систему, соответствующую однородному слою, ее решение дается равенствами (П2.5). Для следующих коэффициентов x(z), y(z) (¿=1,2,3,...) получим дифференциальные уравнения вида
Xi(z) =fi(xo, Уо, xi,yi, x2,y2,..., x¡.¡, yi_1),
¿=1,2,... (П2.10)
yi(z) = gr(xo,yo, xi, У1, x2,y2, ..., xi.i,yi.i),
где f, gi - функции, которые находятся в процессе выполнения описанных преобразований. Так как x0(z), y0(z) удовлетворяют краевым условиям (П2.3), то из (П2.9) следует, что функции xi(z), yi(z) (¿=1,2,3...) удовлетворяют нулевым краевым условиям
x,(Ci) = уЮ = 0, i = 1, 2, 3,... . (П2.11)
Из уравнений (П2.10) следует, что xi(z), yi(z), i=1,2,3... можно последовательно найти простым двукратным интегрированием правых частей (в которых предыдущие коэффициенты ранее найдены) и определением получающихся двух постоянных из условий (П2.11). Для x1(z), y1(z) получим
z Ь
x(z) = (1+b2+c2)/ho [ JMx(s)ds - (z- Q /(¿2- Q) \ Mx(s)ds\
Ь Q (П2.12)
z £
y1(z) = (1+b2+c2)/ho [ \ My(s)ds - (z- Q /(C2- C1) iMy(s)ds],
£
где
b = (&- £) /(&- C1), c = (Л2- m) /(&- C1),
z z
Mx(z) = \(dm /dx)ds, My(z) = \(dm /dy)ds.
<1 <1
В этих равенствах аргументом функции m и ее производных является (xo(s), yo(s), s). Заметим, что в силу принципа стационарности нахождение времени t с точностью до s' требует нахождения траектории (функций x(z), y(z)) с точностью до s'-1. Отсюда, в частности, следует, что для нахождения времени t с точностью до o(s2) достаточно найти x(z), y(z) для '=0,1.
Время t(Q1,Q) пробега волна между точками Q1(^1, tj1, Q) и Q2(%2, т]2, £2) дается равенством
£
t(Q1,Q) = i [1 + x'2(z) + y'2(z)]2 n(x(z), y(z), z) dz.
(1
Подставляя в данное равенство x(z), y(z), найденные по формулам (П2.9) и разлагая его в ряд по степеням s, получим явное приближенное представление времени t в виде
t(Q1, Q2) = to(Q1, Q2) + t1(Q1, Qi)s+ t2(Q1, Q2)i + ... . С точностью до o(s2) оно имеет вид
£ £
t(Q1, Q2) = noD + sD/(C2- C1) imdz + e2{[D/(C2-C1)] i [(¿m /^ + (dm /cy)yd]dz +
Ь Ь (П2.13)
42 42
+ [(£- £)/D] i (bx1,+cy1')mdz + (1/2)и^-i/if/D3 i [x1' 2+y/ 2 + (by/- cx^dz}, ¿2 42
где D = ((£- 6)2 + (Л2- Л1)2 + (z2-zi)2)112 - расстояние между точками Q1, Q2.
Формула (П2.13) дает явное приближенное представление времени t пробега волны в неоднородном слое из точки Q1 в точку Q2 для модели с медленностью (П2.7). Для того, чтобы из данного равенства получить представление для времени пробега волны в среде с медленностью n(x,y,z), достаточно в нем заменить sm(x,y,z) на n(x,y,z)- no. Численные расчеты показали высокую точность данной формулы во многих практически важных случаях. При изменении скорости на 10^15% в области, содержащей точки Q1, Q2, погрешность представления не превосходит 0,1^0,15%.
Заметим, что учет членов порядка s соответствует нахождению времени t вдоль прямой, соединяющей точки Q1, Q2, то есть нахождению времени вдоль траектории, полученной в среде нулевого приближения - однородной среде. В самом деле,
42 42
n,D+ sD/(z2-z1) fm dz = f (1+x'o2 + y'o2)112(no+ sm)dz. i2 i2
Таким образом, первые два слагаемые в (П2.13) дают время пробега волны без учета кривизны луча в неоднородном слое, а третье слагаемое (член порядка s2) позволяет учесть влияние кривизны луча на время пробега волны.
В формулу (П6.7) входят интегралы от медленности и ее производных, при численных же расчетах для конкретных моделей желательно иметь явные формулы, не содержащие интегралы. Если медленность описывается полиномом (или сплайном), зависящим от переменных x, y, z, то интегралы, входящие в равенство (П6.7), берутся в явном виде.
Литература
Goldin S.V. Seismic travel time inversion. Tulsa, 364p., 1986.
Sattleger I.W. A method of computing true interval velocities from expanding spread data in the case of arbitrary long spread and arbitrary dipping plane interfaces. Geophysical Prospecting, v.13, №2, p.306-318, 1975.
Бляс Э.А. Временные поля отраженных волн в трехмерных слоистых со слабо криволинейными границами раздела и латерально-неоднородными слоями. Математические проблемы интерпретации данных сейсморазведки. Новосибирск, Наука, с.98-128, 1988. Бляс Э.А. Кинематика отраженных волн в трехмерных поперечно-изотропных горизонтально-неоднородных средах. Методы расчета и интерпретации сейсмических волновых полей. Новосибирск, Наука, с.148-188, 1991.
Бляс Э.А. Приближенный способ нахождения траекторий лучей в трехмерных слоисто-однородных средах. Геология и геофизика, №12, с.80-87, 1985.
Бляс Э.А., Левит А.Н. Исследование сходимости итерационного алгоритма определения пластовой модели по данным метода ОГТ. Геология и геофизика, №3, с.112-120, 1989.
Глоговский В.М. Анализ методов решения обратной кинематической задачи MOB в неоднородных средах. М., Труды XXXМеждународного геофизического симпозиума, т.2, с.106-116, 1985.
Глоговский В.М., Гогоненков Г.Н. Сходимость итеративного метода определения пластовых скоростей по сейсмическим данным. М., Недра, Прикладная геофизика, вып.92, с.65-78, 1978.
Гольдин C.B. Двумерная кинематическая интерпретация сейсмограмм в слоистых средах. Новосибирск, Наука, 208с., 1993.
Гольдин C.B. Интерпретация данных сейсмического метода отраженных волн. М., Недра, 344с., 1979.
Гольцман Ф.М. Статистические модели интерпретации. М., Недра, 328с., 1971.
Карасик В.М. Изучение скоростей сейсмических волн комплексом методов. М., Недра, 224с., 1993.
Колмогоров А.Н., Фомин C.B. Элементы теории функций и функционального анализа. М., Наука, 324с., 1968.
Маловичко A.A. Кинематическая интерпретация данных цифровой сейсморазведки в условиях вертикально-неоднородных сред. Свердловск, УрО АН СССР, 270с., 1990.
Мешбей В.И. Методика многократных перекрытий в сейсморазведке. М., Недра, 264с., 1985.
Пузырев H.H. Временные поля отраженных волн и метод эффективных параметров. Новосибирск, Наука, 294с., 1979.
Урупов А.К., Левин А.Н. Определение и интерпретация скоростей. М., Недра, 290с., 1985.
Черняк B.C., Гриценко С.А. Интерпретация эффективных параметров ОГТ для пространственной системы однородных слоев с криволинейными границами. Геология и геофизика, №12, с.112-119, 1979.
Яновская Т.Б., Порохова Л.Н. Обратные задачи геофизики. Л., Изд-во ЛГУ, 212с., 1983.
