Исходные сейсмические поля в плосковолновом представлении Текст научной статьи по специальности «Физика»

УДК 550.834.05

ИСХОДНЫЕ СЕЙСМИЧЕСКИЕ ПОЛЯ В ПЛОСКОВОЛНОВОМ ПРЕДСТАВЛЕНИИ

© 2011 г. Г.В. Калайдина, Ю.Д. Борисенко

Кубанский государственный университет, Kuban State University,

ул. Ставропольская, 149, г. Краснодар, 355040, Stavropolskaya St., 149, Krasnodar, 355040,

[email protected] [email protected]

Рассматривается моделирование отражённых продольных PP-, обменных PS- и поперечных SS-волн в докритическом диапазоне в t-p-представлении. В t-p-области волновое поле представляет собой совокупность плоских волн, и в приближении плоскослоистой среды годографы отражённых волн являются эллипсами. Модельные сейсмограммы используются при решении обратных задач сейсморазведки направленными методами Монте-Карло, к которым относятся алгоритм имитации кристаллизации и генетические алгоритмы. Коэффициенты отражения различных типов волн находятся из системы уравнений Цёпприт-ца. Искомые T-p-сейсмограммы получаются в результате свёртки импульсной сейсмограммы с импульсом сейсмического источника.

Ключевые слова: T-p-представление, обменные волны, коэффициент отражения, импульсная сейсмограмма, свёртка, импульс сейсмического источника.

Suggested paper is devoted to the questions of modeling of reflected longitudinal PP-, converted PS- and transverse SS-waves in subcrit-ical band in the т-p-representation. In the т-p-domain the wavefield represents the collection of plane waves and in the flat-layered medium approximation the reflection time-distance curves are the ellipses. Model seismograms apply in the decision of inverse problems of seismics by means of directed Monte-Carlo methods such as simulated annealing and genetic algorithms. The reflectivity factors of different wave modes are found from Zoeppritz system of equations. The required т-p-seismograms are derived as a result of the convolution of impulse seismogram with seismic source pulse.

Keywords: т-p-representation, converted waves, reflectivity factor, impulse seismogram, convolution, seismic source pulse.

Моделирование сейсмических волновых полей как представляет самостоятельный интерес, так и необходимо при решении обратных задач сейсморазведки, в которых используется сравнение зарегистрированных и синтетических сейсмограмм. К таковым относятся направленные методы Монте-Карло (алгоритм имитации кристаллизации [1, 2] и генетические алгоритмы [3]), в которых прямая задача сейсморазведки решается многократно. В вышеупомянутых работах сопоставление зарегистрированных и синтетических сейсмограмм производится в т-р-области. Преобразование Радона - математический аппарат, который в последнее время стал общепринятым при обработке и анализе сейсмических данных. Обычно он более известен как т-р-преобразование - способ наклонного суммирования или метод скоростного накапливания. т-р-преобразование имеет более прозрачный геометрический и физический смысл, чем традиционное преобразование Фурье. Вместо разложения волнового поля на гармонические компоненты т-р-преобра-зование разделяет волновое поле на плосковолновые компоненты с различными временами прихода при нулевом удалении.

В т-р-преобразовании все трассы /-х-сейсмограмм-мы (/ - двойное время пробега; х - удаление источник -приемник) суммируются вдоль заданного линейного направления (наклона), образуя на выходе одну трассу, соответствующую этому значению р (р = А / дх -наклон). Каждый отсчёт выходной трассы на времени т равен, таким образом, сумме отсчетов вдоль траектории / = т + рх. Заметим, что координатами полученного пространства являются лучевой параметр р (параметр Снеллиуса, величина обратная горизонтальной фазовой скорости) и время т на нулевом удалении (время задержки).

Выбор т-р-представления сейсмической информации относительно традиционной /-х-области определяется его возможностями:

- привести годографы различных типов волн к единой эллиптической форме;

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

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

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

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

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

Рассмотрим падение Р-волны на границу раздела двух твердых сред, следуя обозначениям работы [4]. Пусть падает в среде 1 продольная волна на границу раздела сред 1 и 2 с углом падения вх (рис. 1).

i i z /В,

А0 \ V д/ / <h/ /А,

1 X

2 О !/ \ X Хвг

1 d 2Ф, 1 d2 у.

ДФ.----^ = 0, Azi -^^тт- = 0, где j = 1, 2;

- в среде 1:

- в среде 2:

ф = А0в,шд0 + А1в,шд1\ Zi = B^', Ф2 = A2e'aí2 ; Z2 = B2e^2,

a2 dt2

ß2 dt2

л д2 д2 д2

Д = —+--+--- трехмерный оператор Лапласа.

дх ду дг

Временной множитель ехр(—1ю^ в выражениях Ф], X (I = 1, 2) опущен, так как он сокращается в

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

sinöj sinöj sin02 sin$j sin$2

= P, (2)

Рис. 1. Отражённые и преломленные Р- и Б-волны при падении Р-волны на границу раздела двух твёрдых сред

Для двумерного волнового движения в плоскости хОг можно задать скалярный Ф и векторный х потенциалы смещения ^ = УФ + V х х, х = — X • . где 1, к - единичные орты, соответствующие осям х, у, г (ось у направлена на рис. 1 от нас); д д д

У = — 1 +--. +--к - оператор Гамильтона; УФ -

дх ду дг

градиент скалярного потенциала; Ухх - ротор векторного потенциала; ^ = и • 1 + w • к - двумерный вектор смещения. Тогда можно записать

дФ дх дФ дх и =-+ —, w =---- .

дх дг дг дх

Введем обозначения: ^, вх - углы отражения Б-и Р-волны в среде 1; 32, 02 - углы преломления Б-и Р-волны в среде 2; р., а., Д. - плотность, скорость продольных и скорость поперечных волн в |-й среде (I = 1, 2).

Исходя из введенных потенциалов смещения, запишем скалярный и векторный потенциалы в явном виде:

аг аг а2 Д ß

где p = const - лучевой параметр (горизонтальная медленность, параметр Снеллиуса).

С помощью (2) можно упростить выражения (1) для потенциальных функций: g0 = p(x - z ctgöj),

g = p(x + z ctg^), g/= p(x + z ctg ), g2 = p(x - z ctg02),

g2 = p( x - z ctg^).

Запишем граничные условия:

1) непрерывность нормальных смещений на гра-

дФ ду

нице раздела w |j = w |2 при z = 0 или (---) ^ =

dz dx

дФ ду

= (---) |2 при z = 0, откуда получим

dz dx

(-A0 + Ai) ctg в - Bx = A2 ctg в2 - B2;

2) непрерывность сдвиговых смещений дает

.дФ ду дФ ду . U |1 = u|2 пРи Z = 0 или (— + = (— + ^Г)|2 пРи

дк дz дx дz z = 0, откуда найдем A0 + A1 + B1 ctg 31 = A2 - B2 ctg &2;

3) непрерывность нормальных напряжений дает ^ |1 = |2 при z = 0, где ^ = №■\ + 2^ezz, V-\ -дивергенция вектора смещения, которая представляет собой дилатацию (относительное изменение объема);

е22 = — - нормальная деформация; Я и ¡л - парамет-

dz

ры Ламе (л - модуль сдвига).

С учётом вышеприведенных выражений непрерывность нормальных напряжений, принимая во внимание V - (V х х) = 0, примет вид

где A0, A j, Bj - амплитуды потенциалов смещения; ю = 2nv - угловая частота (v - частота волны); i -мнимая единица (i2 = — 1) и g0 =(x sin^ — z cos^V«!, g1 = (x sin^l + z cos^1]/a1,

g = fx sin ^ + z cos 5J/ Д, (1)

g2 =(x sin02 — z cos )/cn2, g2 =(x sin^2 — z cos^2)/ p2.

Потенциалы являются решениями следующих волновых уравнений:

д2Ф д2у 2ДФ + 2И(—ф-—^) dz dxdz

Ii =

—2Ф d2z 2ДФ + 2И{ )

dz dxdz

при z = 0, откуда получим

Л (А + Ai )(1 + ctg2 01) + 2М [(А, + Ai) ctg2 01 - B ctg^i ] =

+

= L2 A2 (1 + ctg2 в2) + 2ju2 (A2 ctg2 в2 + B2 ctg &2), или с учётом выражений для скоростей продольных и

[К

и обобщенного закона Снеллиуса (2), получим

поперечных волн a. =,

\ рИ и ß j =J ^ (/=1, 2)

I2

M(ctg4 - 1)(A0 + A)-2M51ctg51 = = ц2 (ctg2 32 - 1)A2 + 2^2 B2 ctg 32; 4) непрерывность тангенциальных напряжений

CTxz I 1 = CTxz I2 ПРИ Z = 0 где CTxz = 2H^xz, ¿xz =

1 ,du dwч ,

= — (--1--) - сдвиговая деформация.

2 dz dx

С учётом вышеприведенных соотношений получим

d2Ф , d2Z d2/

М(2-+ --f-

dx5z 5z 5x

Ii =

M(2-H--y---\

5x5z 5z dx

откуда найдем ^ [2(-A0 + A1) ctg в1 + B (ctg2 3X -1)] =

= ¡л2 [- 2A2 ctg в2 + B2 (ctg2 32 -1)].

Таким образом, мы получим систему уравнений Кнотта в матричной форме для определения амплитуд потенциалов смещения A, B, A и B2:

ctg 01 1

ft(ctg2 %-1)

- 1

ctg%

-2ft ctg %

ctg 02 -

1

Ctg%2

А

"^2(ctg %2 - 1) - 2ftCtg%2

2ft ctg01 ft(ctg2 %-1)

1- A - A

2ft ctg 0

Л

- ft (ctg 2 %2 -1)

f A) B

A

V B2

(3)

-я(С1Е2 3,-1)• Л0 2 ц Ао ,

После нахождения амплитуд потенциалов смещения можно определить коэффициенты отражения

Грр = А1 / А0 , Гр* = В1 / А0 и прохождения ¡рр = А2 / Ао ,

/ = В /А0 потенциалов.

Чтобы перейти от коэффициентов отражения и прохождения потенциалов к чаще употребляемым коэффициентам отражения и прохождения смещений [5], в уравнениях Кнотта (3) нужно перейти от амплитуд потенциалов смещения к амплитудам смещения:

Ао = а и о, А, = ^П;, В; ; = 1, 2, (4)

ю ю ю

где и ,и - амплитуды смещения в направлении распространения волны; V. - амплитуды смещения по

нормали к направлению распространения волны.

Следует заметить, что точно так же, как в случае с амплитудами потенциалов смещений, амплитуды смещения П, иJ и V. не дают непосредственно амплитуды величин и и V Подставляя (4) в (3), приходим к системе уравнений Цёппритца в матричной форме:

' cos 0 - sin % cos 0 sin % ^ U)

sin01 cos %1 - sin02 cos%2 V

- cos 2% y1 sin 2% m cos 2% my sin2% U 2

/12™201 y cos 2% my22 sin 20 - my cos 2%j V V2 J

cos ei-UQ \ - sin0 - U0 cos 2% - U0 y2 sin 201 -Uo J

где введены обозначения m = — = -Pl2Xl

(5)

Pa1

П =

ß

ß

y2 = —, Z, = p.a. - акустические жесткости (/=1, 2).

Вводя коэффициенты отражения R = Uj IU0, = VIU0 и прохождения

Tpp = U 2IU 0

Tps = VIU0 смещений, система (5) примет вид

( rncfl - ein .Q cos 0

- sin 02 m cos 2% yl2 sin 20 y cos 2% my22 sin 20

cos 0l - sin % sin 0 cos % - cos 2% y sin 2%

sin % cos % my2 sin 2% - my cos 2% j

fR ^

pp

R„.

T

V ps J

cos01

(6)

- sin0j cos 23x

Систему линейных алгебраических уравнений (6) можно решить по правилу Крамера [6]: Rpp=Dpp/D - коэффициент отражения смещений для продольных PP-волн; Rps=Dps/D - коэффициент отражения смещений для обменных PS-волн, где D, Dpp и Dps - определители четвертого порядка, которые соответственно равны

D=

cos0 sin01 - cos 2%

- sin% cos% y1 sin 2%

cos02 - sin02 m cos 2%

y sin20 y cos2% my2 sin20

dpp=

Dps=

cos0 - sin0 cos 2%

- sin% cos% y sin 2%

cos02 - sin02 m cos 2%

y1 sin201 у cos 2% my2 sin202

cos01 cos01 cos02

sin0 - sin0 - sin0

cos 2% cos 2% m cos 2%

y2 sin 20 y2 sin 20 m yl sin 20

sin% cos%2 my sin2% - my cos 2%

sin% cos% my2 sin 2% -my2 cos 2%

sin%2 cos%

шу2

- ту2 cos252

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

Получим выражение пластовой скорости для обменных Р5"-волн. Пусть мощность однородного пласта к, в точке А находится источник продольных Р-волн с пластовой скоростью а (рис. 2).

Рис. 2. К выводу пластовой скорости обменных Р5-волн

В точке В находится приемник поперечных Б-волн, пластовая скорость которых Д В точке О происходит обмен волн: падает продольная Р-волна (луч АО), отражается поперечная Б-волна (луч ОВ).

Здесь введены обозначения: в - угол падения продольной Р-волны; 3 - угол отражения поперечной Б-волны; у = Д/а - отношение скорости поперечной

волны к скорости продольной волны; ур8 - пластовая скорость обменной РБ-волны.

2

а

a

Скорость vps находится как

AO + OB

v„, = -

tp+ts

(7)

где

AO = ■

cos в

OB = . =AO, s = °B (8)

cos 3

a

ß

Подставляя (8) в (7), получим cos в + cos 3

-ja.

(9)

cose + jcos3'

Выражая углы в и 3 через лучевой параметр p из

_ _ sine sin3 обобщенного закона Снеллиуса -=-= p,

соотношение (9) примет вид

yjl - p2a2 +

Г-

2 2 2 j p a

JT-

p2a2 + jj 1 - j2p2a2

a

ja .

ß

(10)

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

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

Аналогично падению Р-волны рассмотрим падение БК-волны (далее обозначаем Б-волна) на границу раздела двух твердых тел под углом 31. В среде 1 под

углом 31 отражается Б-волна, а под углом в1 - Р-вол-на. В среде 2 под углом 32 преломляется Б-волна, под углом в2 - Р-волна. Из граничных условий получим систему линейных алгебраических уравнений Цёп-притца в матричной форме:

(

cose sine

- sin 3 cos 3

cos e - sine

- cos 23 j sin23 m cos 23 j2 sin2e j cos 23 mj22 sin2e i sin3T ^

cos3 j sin 23

sin3

cos 3

mj2 sin23

- mj2 cos 23 y

-jT cos 23

(R Л

sp

R.

Tp

T

\ ps y

(11)

где Rp

U

R. = V-

- коэффициенты отражения сме-

cose sin3

sine cos3

cos23T jTsin23T

jT2sin2eT -jTcos23T

cose - sine

sin32

cos3

m cos 23 mj2 sin 23 mjl sin2e - mj2 cos 23

Рассмотрим плоскослоистую модель среды с числом слоев п. Каждый к-й слой характеризуется плотностью рк, скоростью продольных волн ак, отношением ук = и мощностью hk (к = 1, 2, ..п).

ак

Индексу к = п+1 соответствует подстилающее полупространство с параметрами рк+1, ак+1, Д+1. Для каждой границы раздела в выражениях для коэффициентов отражения Rpp, Rps и Rss углы 31, 32, 0Х и 02 можно выразить через лучевой параметр p с помощью обобщенного закона Снеллиуса (3). Для каждого к-го слоя находим времена задержки тк из уравнения эллипса. Для продольных РР-волн тк , 2ht

где ток = —- - время задержки при р = 0, рак < 1,

ак

так как мы рассматриваем докритический диапазон.

Для обменных волн тк = ток yjT-

2 2

p V psk , где

2h

= ■

V | n =

psk 'p=°

V | n

psk 'p=°

2

psk

находится по формуле (10), а

1 + 7k

jkak, pVps < 1, так как pak < T

щений обменных БР- и поперечных ББ-волн; Т = ,

* V

Т^ = — - коэффициенты прохождения смещений обменных БР- и поперечных ББ-волн.

Из системы (17) найдем коэффициент отражения Я = П , / Д где

Для поперечных SS-волн тк =T0k-J 1 -p2ßk2 , где 2h

ток = —- - время задержки при р = 0, pß < 1, так

ßk

как pak < 1.

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

В качестве сейсмического импульса источника возьмем импульс Н.Н. Пузырева [7]: w(t) = а0 exp(-82i2) х х sin(2^ f0t + ф0), где а0 - амплитудный множитель

(безразмерный параметр); 8- параметр затухания, 1/c; f0 - частота импульса, Гц; t - время, с; <р0 - начальная фаза, рад. Этот импульс имеет «колокольную» огибающую и непрерывен вместе со своими производными. Выбором параметров 8 и ^0 можно получить затухающие колебания с разным характером изменения огибающей. Производя свертку сейсмического импульса с импульсными сейсмограммами, получим г-р-сейсмограммы для различных типов волн.

Создан программный комплекс MDTPWAVE для моделирования г-р-сейсмограмм продольных PP-, обменных PS- и поперечных SS-волн.

На рис. 3 представлены модельные г-р-сейсмограммы для продольных PP-волн (рис. 3а), обменных PS-волн (рис. 3б) и поперечных SS-волн (рис. 3в) шес-тислойной (n = 6) плоскослоистой модели среды, пет-рофизические характеристики которой приведены в таблице.

h

Vps =

Ds=

Р, мс/км Р, мс/км Р мс/км

80 160 240 320 400 80 160 240 320 400 80 160 240 320 400

Рис. 3. Модельные т-р-сейсмограммы для различных типов волн: а - продольные РР-волны; б - обменные PS-волны; в - поперечные Ж-волны

Петрофизические характеристики шестислойной модели среды

к ph г/см3 ак, км/с Ук hk, м

1 1,90 2,50 0,30 50

2 2,00 2,65 0,35 60

3 2,10 2,80 0,40 50

4 2,25 2,90 0,45 40

5 2,20 2,75 0,38 50

6 2,30 3,00 0,50 60

7 2,35 3,30 0,55 -

В программном комплексе были заданы также следующие параметры: А0 = 1000 - амплитуда импульсной сейсмограммы; Тъ = 0 мс - начальное время сейсмограммы, Те = 2000 мс - конечное время сейсмограммы; Д/ = 2 мс - шаг дискретизации по времени; /0 = 35 Гц -частота импульса Пузырёва; 8= 20 1/с - параметр затухания; МТ = 1 - число периодов импульса; а0 = 1000 -амплитудный множитель импульса Пузырева; к = 6 -число слоёв модели; ръ = 4 мс/км - начальное значение лучевого параметра р; ре = 404 мс/км - конечное значение лучевого параметра р; Др = 4 мс/км - шаг по луче-

вому параметру p; Xdp = 1000 м - координата сейсмограммы ОСТ (общая срединная точка).

Модельные T-p-сейсмограммы обменных PS- и поперечных SS-волн могут быть использованы в задачах обращения T-p-сейсмограмм в параметры модели среды (подобно обращению T-p-сейсмограмм продольных PP-волн в параметры модели среды [2]).

Литература

1. Sen M.K., Stoffa P.L. Nonlinear one-dimensional seismic waveform inversion using simulated annealing // Geophysics. 1991. Vol. 56, № 10. P. 1624 - 1638.

2. Курочкин А.Г., Борисенко Ю.Д., Калайдина Г.В. Инверсия сейсмической информации в параметры модели среды // Геофизика. 2003. Спец. вып. «Технология сейсморазведки II». С. 44 - 47.

3. Stoffa P.L., Sen M.K. Nonlinear multiparameter optimization using genetic algorithms: Inversion of plane-wave seis-mograms // Geophysics. 1991. Vol. 56, № 11. P. 1794-1810.

4. Шериф Р., Гелдарт Л. Сейсморазведка : в 2 т. Т. 1 : История и получение данных : пер. с англ. М., 1987. 448 с.

5. Рябинкин Л.А. Теория упругих волн : учеб. пособие для вузов. М., 1987. 182 с.

6. Корн Г., Корн Т. Справочник по математике : пер. с англ. М., 1968. 620 с.

7. Боганик Г.И., Гурвич И.И. Сейсморазведка : учебник для вузов. Тверь, 2006. 744 с.

Поступила в редакцию

2 декабря 2010 г.