Расчет волновых полей для анизотропных сред с поглощением на основе отражательного метода

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

Повний опис

Збережено в:
Бібліографічні деталі
Дата:2011
Автори: Роганов, Ю.В., Роганов, В.Ю.
Формат: Стаття
Мова:Russian
Опубліковано: Інститут геофізики ім. С.I. Субботіна НАН України 2011
Назва видання:Геофизический журнал
Онлайн доступ:http://dspace.nbuv.gov.ua/handle/123456789/97099
Теги: Додати тег
Немає тегів, Будьте першим, хто поставить тег для цього запису!
Назва журналу:Digital Library of Periodicals of National Academy of Sciences of Ukraine
Цитувати:Расчет волновых полей для анизотропных сред с поглощением на основе отражательного метода / Ю.В. Роганов, В.Ю. Роганов // Геофизический журнал. — 2011. — Т. 33, № 4. — С. 117-126. — Бібліогр.: 24 назв. — рос.

Репозитарії

Digital Library of Periodicals of National Academy of Sciences of Ukraine
id irk-123456789-97099
record_format dspace
spelling irk-123456789-970992016-03-26T03:02:20Z Расчет волновых полей для анизотропных сред с поглощением на основе отражательного метода Роганов, Ю.В. Роганов, В.Ю. Запропоновано стійкий спосіб обчислення хвильових полів для горизонтально-шаруватих анізотропних середовищ на основі методу відбиття. Для градієнтних середовищ виведено диференціальне рівняння Ріккаті, коефіцієнти якого визначають за методом збурень. A stable method for calculating wave fields in horizontally layered anisotropic media on the basis of the reflectivity method is offered. For gradient media, a differential Riccati equation with coefficients defined by the perturbation method is derived. Предлагается устойчивый способ расчета волновых полей для горизонтально-слоистых анизотропных сред на основе отражательного метода. Для градиентных сред выведено дифференциальное уравнение Риккати, коэффициенты которого определяются по методу возмущений. 2011 Article Расчет волновых полей для анизотропных сред с поглощением на основе отражательного метода / Ю.В. Роганов, В.Ю. Роганов // Геофизический журнал. — 2011. — Т. 33, № 4. — С. 117-126. — Бібліогр.: 24 назв. — рос. 0203-3100 http://dspace.nbuv.gov.ua/handle/123456789/97099 550.834 ru Геофизический журнал Інститут геофізики ім. С.I. Субботіна НАН України
institution Digital Library of Periodicals of National Academy of Sciences of Ukraine
collection DSpace DC
language Russian
description Запропоновано стійкий спосіб обчислення хвильових полів для горизонтально-шаруватих анізотропних середовищ на основі методу відбиття. Для градієнтних середовищ виведено диференціальне рівняння Ріккаті, коефіцієнти якого визначають за методом збурень.
format Article
author Роганов, Ю.В.
Роганов, В.Ю.
spellingShingle Роганов, Ю.В.
Роганов, В.Ю.
Расчет волновых полей для анизотропных сред с поглощением на основе отражательного метода
Геофизический журнал
author_facet Роганов, Ю.В.
Роганов, В.Ю.
author_sort Роганов, Ю.В.
title Расчет волновых полей для анизотропных сред с поглощением на основе отражательного метода
title_short Расчет волновых полей для анизотропных сред с поглощением на основе отражательного метода
title_full Расчет волновых полей для анизотропных сред с поглощением на основе отражательного метода
title_fullStr Расчет волновых полей для анизотропных сред с поглощением на основе отражательного метода
title_full_unstemmed Расчет волновых полей для анизотропных сред с поглощением на основе отражательного метода
title_sort расчет волновых полей для анизотропных сред с поглощением на основе отражательного метода
publisher Інститут геофізики ім. С.I. Субботіна НАН України
publishDate 2011
url http://dspace.nbuv.gov.ua/handle/123456789/97099
citation_txt Расчет волновых полей для анизотропных сред с поглощением на основе отражательного метода / Ю.В. Роганов, В.Ю. Роганов // Геофизический журнал. — 2011. — Т. 33, № 4. — С. 117-126. — Бібліогр.: 24 назв. — рос.
series Геофизический журнал
work_keys_str_mv AT roganovûv rasčetvolnovyhpolejdlâanizotropnyhsredspogloŝeniemnaosnoveotražatelʹnogometoda
AT roganovvû rasčetvolnovyhpolejdlâanizotropnyhsredspogloŝeniemnaosnoveotražatelʹnogometoda
first_indexed 2025-07-07T04:28:01Z
last_indexed 2025-07-07T04:28:01Z
_version_ 1836960954957430784
fulltext РАСЧЕТ ВОЛНОВЫХ ПОЛЕЙ ДЛЯ АНИЗОТРОПНЫХ СРЕД С ПОГЛОЩЕНИЕМ НА ОСНОВЕ... Геофизический журнал № 4, Т. 33, 2011 117 УДК 550.834 Расчет волновых полей для анизотропных сред с поглощением на основе отражательного метода © Ю. В. Роганов, В.Ю. Роганов, 2010 Украинский государственный геологоразведочный институт, Киев, Украина Поступила 15 февраля 2010 г. Представлено членом редколлегии Ю. К. Тяпкиным Запропоновано стійкий спосіб обчислення хвильових полів для горизонтально-шаруватих анізотропних середовищ на основі методу відбиття. Для градієнтних середовищ виведено диференціальне рівняння Ріккаті, коефіцієнти якого визначають за методом збурень. A stable method for calculating wave fields in horizontally layered anisotropic media on the basis of the reflectivity method is offered. For gradient media, a differential Riccati equation with coefficients defined by the perturbation method is derived. Введение. Моделирование волновых полей является важным инструментом при построе- нии изображений, интерпретации, инверсии сейсмических данных, проектировании систем наблюдений. Благодаря развивающейся ком- пьютерной технологии становятся доступными различные методы расчета волновых полей для сложных геологических сред, основанные на уравнениях динамической теории упругости. Однако синтез волнового поля на основе пло- ских гармоник, рассчитываемых матричным способом, по-прежнему широко использует- ся благодаря его уникальным свойствам. Этот метод позволяет учитывать эффекты анизо- тропии и частотно-зависимого поглощения, наличие свободной поверхности, задавать произвольные типы источников, получать ча- стичные волновые поля, содержащие заданные типы волн. Однако существуют препятствия при непосредственной компьютерной реализа- ции этого метода, вызванные необходимостью вычисления миноров плохообусловленных матриц. При этом требуется находить суммы возрастающих и убывающих экспонент, что приводит к потере точности вычислений [Ро- ганов, 2009]. Были предложены разные методы повы- шения точности вычислений, такие как ис- пользование миноров вместо пропагаторов [Dunkin, 1965], сведение задачи к пятимерному формализму [Молотков, 1984], рекурсивные формулы расчета миноров с учетом различ- ного поведения экспонент в зависимости от горизонтальной медленности [Abo-Zena, 1979]. Наиболее удобным и широко используемым методом стабилизации вычислений является отражательный метод, предложенный К. Фук- сом и Г. Мюллером [Fuchs, Muller, 1971]. В дальнейшем этот метод интенсивно развивал- ся Б. Л. Кеннетом [Kennett, 1975; 1979; 1983], Р. Кайндом [Kind, 1976], Г. Фраером [Fryer, 1981], С. Маликом и Н. Фразером [Mallick, Frazer; 1987; 1988] и др. В начальной версии метод был создан для пачки упругих слоев, расположенных между упругими полупространствами с источником и приемником, находящимися над областью рас- сеяния. Р. Степхен [Stephen, 1977] и Б. Л. Кен- нет [Kennett, 1975] модифицировали его для различных типов источников и приемников, расположенных внутри слоев под свободной поверхностью. Учет поглощения в теории был выполнен в работах [Kennett, 1975, Sipkin et. al; 1978, O’Neil, Hill, 1979]. Отражательный метод был применен для расчета волновых полей в анизотропных [Fryer, Frazer, 1984] и градиентных средах [Mallick, Frazer, 1987], для различных систем наблюдений, включая ВСП [Mallick, Frazer, 1987; 1988]. Суть метода состоит в последовательном добавлении слоев к пачке в процессе итера- тивного расчета матриц рассеяния. Учитывая, что возрастающие компоненты не содержатся в матрицах рассеяния и не возникают в итера- ционной процедуре, этот способ расчета коэф- фициентов отражения и преломления является устойчивым для всех частот и горизонтальных медленностей. Ю. В. РОГАНОВ, В.Ю. РОГАНОВ 118 Геофизический журнал № 4, Т. 33, 2011 Поскольку каждая частота и горизонталь- ная медленность обрабатываются независи- мо, алгоритм расчета хорошо распараллели- вается, позволяя получить трехкомпонентное 3D-волновое поле без интенсивных временных затрат. В настоящей статье отражательный рекур- сивный метод объединяется с методом воз- мущений для получения дифференциального уравнения Риккати, которому удовлетворяет матрица рассеяния в градиентной анизотроп- ной среде, а также изучаются свойства раз- личных точечных источников в анизотропных средах. Теория. Рассмотрим точечный источник, находящийся в начале системы координат x=0, состоящий из однонаправленной силы f=(fi), пары сил ( )A A k ikM=M и двойной пары сил ( )S S k ikM=M . Выражения A kM , S kM и Mk определяют k-е столбцы соответствующих матриц. Обозначим A S k k k= +M M M полный момент сил. Известно, что матрица пары сил A ikM является кососимметрической, а матри- ца двойной пары сил S ikM — симметриче- ской, т. е. выполняются равенства A A ik kiM M  ,S S ik kiM M= . Распространение упругих волн, возбужден- ных этим источником в анизотропной упру- гой среде, описывается дифференциальными уравнениями ( ) ( ) k Ai ik i ik x k u f M t x x x , (1) ( ), S pim im im pk k u M t x t x , (2) где u=(ui) — вектор скоростей смещений, τk=(τik) — тензор напряжений, λij,mn — тензор упругих постоянных, ρ — плотность. По повторяющим- ся индексам k, p=1,2,3 проведено суммирование. Предполагая, что свойства среды не изме- няются вдоль латерали, выполним 3D-преобра- зование Фурье дифференциальных уравнений (1), (2) по переменным (x1,x2,t): ( ) ( ) ( )1 1 2 2 1 2 1 2 1 2, , , , j t x p x pg p p g x x t e dx dx dt . Обратное Фурье-преобразование описыва- ется соотношением ( ) ( ) ( )1 1 2 22 1 2 1 2 1 2, , , , j x p x p tg x x t g p p e dp dp dt и позволяет синтезировать волновое поле по плоским гармоникам. Учитывая, что в спек- тральной области производным соответству- ют умножения на множители 1 1 j p x , 2 2 j p x , j t , получим соотношения ( )3 3 3 n nj j p x x u f ( ) ( )3 3 3 3 A A n n x j p x x M M , (3) m n mnj j p C u ( )3 3 3 S m mj x x uC M . (4) Здесь и в дальнейшем предполагается, что по повторяющемуся индексу n=1,2 выполне- но суммирование, а тензор упругих посто- янных представляется в виде набора матриц (λij,mn)=Cim[j,n]. Систему уравнений (3)—(4) можно предста- вить в виде матричного дифференциального неоднородного уравнения, если выполнить следующие алгебраические преобразования. Из уравнения (4) при m=3 найти значение ∂u/∂x3 и подставить его в уравнение (4) при m=1, 2. Затем полученные значения τ1 и τ2 из уравнения (4) подставить в уравнение (3) и из него найти ∂ 3/∂x3. В результате получается уравнение 3 d j dx v Fv g , (5) где 3 = u v — шестимерный вектор скоростей смещений-напряжений, 1 33 T A C F B A , ( )1 33 1 31 2 32p p= +A C C C , ( )1 3 33 3 , 1,2 m n m n mn m n p p = B C C C C I , ( )3 3 000 n S n j p j xg F Mf M ( ) 3 3 3 0 xA x M . (6) I — единичная 3×3-матрица. Уравнение (5) имеет решение РАСЧЕТ ВОЛНОВЫХ ПОЛЕЙ ДЛЯ АНИЗОТРОПНЫХ СРЕД С ПОГЛОЩЕНИЕМ НА ОСНОВЕ... Геофизический журнал № 4, Т. 33, 2011 119 ( ) ( ) ( ) 3 3 3 0 exp exp x x j x j dv F F g . (7) Подставляя в (7) значение g из равенства (6), находим, что источник создает разрыв скоро- стей смещений и напряжений Δv в плоскости x3=0: 3 0 00 n n j p jv F M Mf . (8) Аналогичное соотношение получено в рабо- тах [Kennett, 1983; Fryer, Frazer, 1984]. Отметим, что при выводе соотношения (8) учитывалось, что A S n n n+ =M M M , а также ( ) ( ) 3 3 30 0 0 exp x A Aj d jF F M M . Собственными числами матрицы F являют- ся вертикальные медленности ξα всех типов плоских волн, распространяющихся вверх и вниз с горизонтальными медленностями (p1,p2). Собственными векторами этой матри- цы являются векторы скоростей смещений- напряжений vα=(aα,τα)T, которые нормируем так, чтобы 1=a . Обозначая символами d и u вниз и вверх распространяющиеся волны, а qP, S1 и S2 — их типы, объединим векторы vα и вертикальные медленности ξα в матри- цы размера 6×6: ( )1 2 1 2, , , , ,d d d u u u qP S S qP S S=E v v v v v v   , ( )1 2 1 2diag , , , , ,d d d u u u qP S S qP S SL . Пусть 0 0 = I J I . Поскольку матрицы B и 1 33C симметричны, то JFJ=FT. Поэтому собственными левыми вектор-строками матрицы F являются векторы ( ), 2 2 T TT T T= = av J w a a , где α — типы всех волн, распространяющихся вверх и вниз. Эти век- торы являются строками матрицы E–1. Разрыв скоростей смещений и напряжений Δv эквива- лентен разрывам амплитуд нисходящих и вос- ходящих плоских волн Δs=(δd–δu) T=E–1Δv. Вос- ходящая волна типа α имеет разрыв амплитуды 3 2 2 2 T T T n nT T Td j p j a f a M a M w v a a a . Аналогичная формула со знаком минус справедлива для нисходящих волн. При вы- воде формулы использовалось соотношение wαF=ξαwα. Обозначим nα единичный вектор, направ- ленный вдоль распространения волны типа α; να — соответствующую фазовую скорость. По- скольку ( )1 2 1, , Tp p v n , то 2 2 T T T Td j v a f a Mn a a . (9) Для изотропной среды 2T va и, следо- вательно, 2 32 2 T T d j v v a f a Mn . (10) Формула (10) обычно выводится с использова- нием теоремы Ламе и свойств функции Грина, справедливых только для изотропной среды [Аки, Ричардс, 1983]. Отметим, что формулы (9), (10) определяют характеристики направ- ленности для точечных источников разных типов в однородной среде и являются множи- телями при синтезе волнового поля. Рассмотрим среду, состоящую из полупро- странства i=n+1, над которым расположены n однородных слоев с мощностями hi (i=1,…,n) и плоскими горизонтальными границами раз- дела. Сверху над слоями расположено либо полупространство i=0, либо свободная поверх- ность (рис. 1). Обозначим плотности и пара- метры упругости слоев и полупространств со- ответственно ρi и ( ) , i mp nq (i=0,…,n+1). Полупро- странства и слои могут быть анизотропными с наличием поглощения. Предполагается, что слои находятся в жестком контакте. Будем ис- пользовать декартовую систему координат с осью OX3, перпендикулярной границам раздела и направленной вниз. Обозначим zi=h1+…+hi глубину залегания нижней границы i-го слоя. Рис. 1. Краевые условия и обозначения для пачки слоев между упругими полупространствами (а) и свободной по- верхностью и упругим полупространством (б). Ю. В. РОГАНОВ, В.Ю. РОГАНОВ 120 Геофизический журнал № 4, Т. 33, 2011 Если между глубинами zk и zm нет раз- рывных источников, векторы скоростей смещений-напряжений v(z) являются непре- рывными функциями от z на интервале zm ≥ z  ≥ zk и связаны соотношением v(zm)=P(zm,zk)v(zk), где матрица ( ),m kz z =P ( ) ( )1 1exp ...expm m k kj h j h + +F F (11) является пропагатором пачки слоев, Fi — ма- трица i-го слоя, входящая в уравнение (5). Про- пагатор слоя P(zi,zi–1) можно представить в виде ( ) 1 1,i i i i iz z =P E E . (12) Между вектором амплитуд b(z) всех типов волн на глубине z в i-м слое и вектором скоро- стей смещений и напряжений v(z) выполняется соотношение v(z)=Eib(z). Поэтому векторы ам- плитуд всех типов волн на глубинах kz и mz + связаны равенством ( ) ( ) ( ),m m k kz z z z=b Q b , где ( ) ( )1 1, ,m k m m k kz z z z+=Q E P E — амплитудный пропагатор пачки слоев. Поскольку функция b(z) имеет разрывы на глубинах расположения границ zi, необходимо различать точки над и под границей, обозначая их верхними индек- сами – и + соответственно. Используя равенство (12), получаем ( ) 1 1, ...m k m m k k kz z + +=Q G G G , (13) где 1 1i i i+=G E E — матрица i-й границы. Отме- тим, что матрицы Λi — диагональные а Gi — частотно-независимые для среды без погло- щения. Амплитудный пропагатор ( ),m kz zQ связывает между собой амплитуды всех типов волн над k-й под m-й границами. Аналогично матрица рассеяния ( ),m kz zS связывает ампли- туды волн, направленных к пачке слоев zm≥z≥zk, с амплитудами рассеянных этой пачкой волн и направленных во внешнюю сторону. Разбивая векторы b на трехкомпонентные подвекторы, соответствующие нисходящим и восходящим волнам, и представляя матрицы ( ),m kz z=Q Q и ( ),m kz z=S S в блочном виде, получим ( ) ( ) ( ) ( ) 11 12 21 22 d m d k u m u k z z z z = b bQ Q Q Qb b , (14) ( ) ( ) ( ) ( ) 11 12 21 22 d m d k u k u m z z z z = b bS S S Sb b . (15) Заметим, что переход от матрицы Q к ма- трице S состоит в перестановке векторов ( ) ( )u m u kz zb b в равенствах (14)—(15). По- скольку такая перестановка, применяемая дважды, — тождественный оператор, то фор- мулы пересчета Q в S такие же, как формулы пересчета S в Q, просто в них необходимо об- менять местами S и Q: 1 11 11 12 22 21S Q Q Q Q , 1 12 12 22=S Q Q , 1 21 22 21S Q Q , 1 22 22=S Q , (16) 1 11 11 12 22 21Q S S S S , 1 12 12 22=Q S S , 1 21 22 21Q S S , 1 22 22=Q S . (17) По смыслу матрицы S11, S21 содержат все коэффициенты преломления и отражения со- ответственно при инициализации пачки сло- ев сверху. Поэтому их обозначают TD и RD. Аналогично матрицы S12 и S22 содержат ко- эффициенты отражения и преломления при инициализации пачки слоев снизу. Их обозна- чают RU и TU. Основная идея отражательного метода со- стоит в применении формул для расчета матри- цы рассеяния композита двух пачек слоев по матрицам рассеяния каждой из пачек. Вывод этих формул основан на определении ампли- тудных пропагаторов пачек слоев по формулам (17), перемножении найденных пропагаторов и вычислении матрицы рассеяния по полученно- му произведению с использованием формулы (16): ( ) ( ) ( ) ( ) ( ) ( ) , , , . , , , . m k m k m k m k z z z z z z z z z z z z S S S Q Q Q . (18) Операция получения матрицы рассеяния композита двух пачек слоев называется звезд- ным произведением [Ursin, 1983]. Из соотноше- ния (18) следует, что коэффициенты прелом- ления и отражения, содержащиеся в матрице рассеяния, можно вычислить по следующим формулам [Kennett, 1975; Kennett, Kerry, 1979]: ( ) 1 D D U D DT T I R R T , (19) ( ) 1 U U D U D U UR R T R I R R T , (20) ( ) 1 U U D U UT T I R R T , (21) РАСЧЕТ ВОЛНОВЫХ ПОЛЕЙ ДЛЯ АНИЗОТРОПНЫХ СРЕД С ПОГЛОЩЕНИЕМ НА ОСНОВЕ... Геофизический журнал № 4, Т. 33, 2011 121 ( ) 1 D D U D U D DR R T R I R R T . (22) В соотношения (19)—(22) входят матрицы, определяемые равенствами ( ), D U m k D U z z = T R S R T , ( ), D U m D U z z+ = T R S R T , ( ), D U k D U z z = T R S R T . В равенствах для RU и RD можно изменить порядок сомножителей, пользуясь формулами ( ) ( )1 1 U D U U D UR I R R I R R R , ( ) ( )1 1 D U D D U DR I R R I R R R . При выводе формул (19)—(22) не исполь- зуется конкретный вид элементов матриц Q. Поэтому они справедливы для любого разло- жения амплитудного пропагатора в произве- дение. В частности, матрицу рассеяния мож- но вычислить с использованием разложения (13) по формуле ( ),m kz z =S ( ) ( )* *...m mS G S ( ) ( ) ( )1 1* * *k k k+ +S G S S G . Аналогично приме- няется разложение (11) и равенство ( ),m kz z =S = ( )1 1 *m+S E ( ) ( )( , ) *m k kz zS P S E . Экспоненциаль- ная зависимость от частоты присутствует толь- ко в матрицах i=diag( D, U), причем для них TD= D, 1 U U=T , RD=RU=0. Матрицы TD,TU, RD иRU не содержат экспонент с положитель- ными показателями и композиция с ними по формулам (19)—(22) также не изменяет это свойство. Отсюда следует, что в процессе вы- числения коэффициентов матрицы рассеяния не возникают возрастающие экспоненты. Сле- довательно, процесс вычисления коэффици- ентов отражения и прохождения на основе отражательного метода является устойчивым для всех медленностей и частот [Kennett, 1983]. Определим волновое поле u(zr), зарегистри- рованное приемниками в точке z=zr, возбуж- денное точечным источником, находящимся на глубине z=zs в одном из слоев пачки с ис- пользованием только матриц рассеяния. Для пачки слоев, расположенных между упругими полупространствами, выполняются условия излучения ( )0 0d z =a и ( ) 0u nz + =a , а волновое поле имеет разрыв амплитуд s=( d– – u) T на глубине zs. Следовательно, краевая за- дача описывается уравнением ( ) ( )0 0 0 ,n u z z z +Q a ( ) ( ), 0 d d n n s u z z z + ++ = a Q . (23) Для пачки слоев, расположенных под сво- бодной поверхностью, справедливы соотноше- ния ( )3 0 0z+ = и ( ) 0u nz + =a . Поэтому краевая задача сводится к уравнению ( ) ( ) 1 0 1 0 0 ,nz z z+ +Q E u ( ) ( ), 0 d d n n s u z z z + ++ = a Q , (24) где ( ) ( ) ( ) ( ) 0 3 0 1 0 0 d u z z z z + + + + = a E a u , т. е. 21 22 1 11 12 = E E E E E   , и Eij (i,j=1,2) — блоки матрицы собственных вектор ов слоя под свободной поверхностью. Из соотношений (23), (24) следует, что крае- вые задачи для пачки слоев, расположенных между упругими полупространствами, и пачки слоев под свободной поверхностью описыва- ются уравнениями одинакового типа, только в первом случае матрица верхней границы 1 0 1 0=G E E , а во втором случае 1 0 1=G E . Из соотношений (23), (24) определяются значения ( )0u za и ( )d nz +a для задачи с упруги- ми полупространствами либо значения ( )0z +u и ( )d nz +a для задачи со свободной поверхностью. C использованием полученных данных опреде- ляется искомый вектор ( ) ( ) ( )3 r r r z z z = u v . Фор- мулы получаются разными и зависят от взаим- ного расположения источника и приемника. Если приемник находится под источником, выполняется пересчет вверх вектора амплитуд с глубины nz z+= на глубину rz z= : ( ) ( ) ( )1 , 0 d n r r n r z z z z + = a v E Q . (25) Если приемник находится над источником, выполняется пересчет вниз либо вектора ам- плитуд с глубины 0z z= на глубину z   =   zr для пачки слоев между упругими полупростран- ствами: Ю. В. РОГАНОВ, В.Ю. РОГАНОВ 122 Геофизический журнал № 4, Т. 33, 2011 ( ) ( ) ( )0 0 0 ,r r r u z z z z =v E Q a , (26) либо вектора скоростей смещений и напряже- ний для пачки слоев под свободной поверхно- стью 0z z+= : ( ) ( ) ( ) 1 0 1 0 0 ,r r rz z z z+ =v E Q E u . (27) При использовании отражательного метода все пропагаторы Q и матрицы E, входящие в формулы (23)—(27), выражаются через соот- ветствующие матрицы рассеяния по формулам (16), а произведения заменяются соотношения- ми (19)—(22). Обозначим для краткости глубины z0, zs, zr, zn их индексами 0, s, r, n соответственно. После выполнения алгебраических преобразований получаются формулы для вектора скоростей смещений u(zr) волн в пачке слоев между упру- гими полупространствами, если приемник рас- положен выше источника: ( ) ( ) ( ) 10 0 12 11 r r r sr r r u d uzu E E R I R R ( ) ( )10sr ns s ns u d u u d dT I R R R , (28) и ( ) ( ) ( ) 1 12 11 r nr r rs nr rs r d u d dzu E R E I R R T ( ) ( )10 0s ns s u d d u uI R R R , (29) если приемник расположен ниже источника. В формулах (28) и (29) приняты сокращения вида ( ),nr d d n rz z+=R R . При наличии свободной поверхности фор- мулы (28) и (29) остаются в силе, только матри- цы отражений 0s uR и 0r uR следует вычислять соответственно по пропагаторам ( ) 1 0 1,sz zQ E и ( ) 1 0 1,rz zQ E . Обозначим сверху значком «~» матрицы отражения и преломления для среды с продолженным вверх первым слоем, т. е. без свободной поверхности. Эти матрицы вычис- ляются по пропагатору ( )0,z z+Q . Обозначим также ( ) ( )0 1 11 0 12 0 f u z z=R E E . Если воспользо- ваться соотношением (20), можно вывести формулы ( ) 10 0 0 0 0 0 0s s s f s f s u u d u d u uR R T R I R R T , ( ) 10 0 0 0 0 0 0r r r f r f r u u d u d u uR R T R I R R T . Для получения вектора скоростей смеще- ний выполняется обратное преобразование Фурье: ( )1 2, , ,rx x z t =u (30) ( ) ( )( )2 1 1 2 2 1 2exprz j x p x p t dp dp du . Уменьшение искажений, связанных с яв- лением Гиббса, возникающим при интегри- ровании в конечной области медленностей 2 2 2 1 2 maxp p p+ < , выполняется умножением под- ынтегральной функции на окно r(p1,p2), полу- ченное билинейным преобразованием АЧХ фильтра Баттерворта [Каппелини и др., 1983]: 0 2tg 5 f = , 2 2 1 2 max tg 2 p p f p = ; ( ) ( ) 1 2 81 0 1, 1 r p p ff = + , где 1 max min1,2p v= , a νmin — наименьшая скорость в слоях. При расчете волновых полей градиентные зоны необходимо разбивать на тонкие слои. Альтернативным подходом является решение дифференциального уравнения Риккати, ко- торому удовлетворяет матрица рассеяния в градиентной зоне [Abramovici, 1968]. Выведем это уравнение. Пропагатор P(z,z0) как функция глубины z удовлетворяет дифференциальному уравне- нию d dz jP FP , (31) где F=F(z) — матрица слоя. Подставляя в (31) выражение ( ) ( ) ( ) ( )1 0z z z z=P E Q E , получaeм уравнение для амплитудного пропагатора [Ursin, 1983; Kennett, 1983]: 1 dd dz j dz EQ L E Q , (32) где F=ELE–1, L=L(z) — диагональная матри- ца содержащая вертикальные медленности; E=E(z) — матрица, составленная из собствен- ных вектор-столбцов матрицы F(z). Разложим столбы матрицы dE/dz по столб- цам матрицы E, т. е. воспользуемся формулой d dz =E ED , (33) где D=E–1dE/dz — матрица коэффициентов раз- ложения. Подставляя соотношение (33) в (32), получaeм уравнение для амплитудного пропа- гатора Q=Q(z, z1) в виде ( )d dz jQ L D Q . (34) РАСЧЕТ ВОЛНОВЫХ ПОЛЕЙ ДЛЯ АНИЗОТРОПНЫХ СРЕД С ПОГЛОЩЕНИЕМ НА ОСНОВЕ... Геофизический журнал № 4, Т. 33, 2011 123 Пропагатор Q=Q(zn, z) как функция от верх- ней глубины удовлетворяет аналогичному уравнению с переставленными множителями в правой части: ( )d dz jQ Q L D . (35) Значения элементов матрицы D=(dαβ) опре- деляются зависимостью F(z). Чтобы их найти, воспользуемся теорией возмущений, предпо- лагая, что матрица F(z) приводится к диагональ- ному виду подобными преобразованиями. Дифференцируя равенство Fvβ=ξβvβ и умно- жая результат слева на wα, получим ( )w F v w v w v . Поэтому, если α≠β и ξα≠ξβ, то d = = w F v w v . Если ξα=ξβ, то остается свобода в выборе собственных векторов vα и vβ матрицы F. Их можно выбрать так, что dαβ=0, при α≠β. Значе- ние dαα определяется условием нормировки ( ) 1z =a , из которого следует ( )1 Td d a a . Уравнения (34) и (35) являются частными случаями уравнения d dzQ MQ QN , (36) гдеM=j L–D, N=0, либо M=0, N=j L–D. ПустьM=(Mij), N=(Nij) — представление ма- триц в блочном виде и S=(Sij) — матрица рас- сеяния, соответствующая пропагатору Q=(Qij). Покажем, что матрица S удовлетворяет ма- тричному уравнению Риккати: 12 11 21 22 0 0 0 0 d dz M MS S N N 11 12 22 21 0 0 0 0 N N S S S M M . (37) Для этого рассмотрим матрицы 11 1 21 0 = S F S I , 12 2 220 = I S F S , 11 12 3 0 = S S F I , 4 21 22 0 = I F S S , 5 21 22 0 = I F Q Q . Используя равенства (16) и (17), несложно показать, что 1 2 0+ =F F Q , (38) 4 3=QF F , (39) 4 5 =F F I . (40) Подстановкой можно также проверить, что 1 2 5 d dd dz dz dz F FS F Q . (41) Дифференцируя равенство (38) и подстав- ляя в результат правую часть равенства (41), находим 5 2 d d dz dz = S QF F , т. е. 2 4 d d dz dz = S QF F . Умножая равенство (36) слева на F2 и справа на F4, получаем 2 3 1 4d dz = +S F MF F NF . Расписывая M и N в блочном виде, получaeм соотношение (37). Иногда рассматривают матрицу рассеяния R с блочным видом [Ursin, 1983]: 21 22 11 12 11 12 21 22 0 0 = = S S S SI R S S S SI . Уравнение Риккати для матрицы R получа- ется, если умножить соотношение (37) слева и справа на матрицу J и учесть, что J2=I: 21 22 12 11 0 0 0 0 d dz N NR R M M 11 12 22 21 0 0 0 0 N N R R R M M . В таком виде уравнение Риккати пред- ставлено в статье [Ursin, 1983] для изотропной слоистой среды. В статье [Norris, Shuvalov, 2010] выведено уравнение Риккати, которо- му удовлетворяет двухточечная импедансная матрица Z(zm, zk). Блоки этой матрицы можно получить по формулам (16), примененным к прорагатору P(z, z0). На основе импедансной матрицы легко получить матрицу рассеяния, выполняя звездные умножения ( ),m kz z =S ( ) ( ) ( )1 1 * , *m m k kz z+= S E Z S E . Численный пример. Применение отража- тельного метода для формирования волнового поля продемонстрируем на примере однород- ной среды с вертикальной системой трещин, Ю. В. РОГАНОВ, В.Ю. РОГАНОВ 124 Геофизический журнал № 4, Т. 33, 2011 перечных волн VSV=1300 м/с и плотностью ρ=1200 кг/м3. Трещиноватость моделируется в двух вариантах: без поглощения и с поглоще- нием. В последнем случае она задается ком- плексными значениями ослабленностей ΔN и ΔT [Chichinina et al., 2009]. Схема наблюдений состоит из одного источ- ника в начале системы координат и 19 приемни- ков, расположенных вдоль дуги окружности на расстоянии 1 км от источника в вертикальной плоскости, содержащей ось OX (рис. 2). Углы α между направлениями источник-приемник и вертикалью изменяются в пределах от 0 до 90° с шагом 5°. На рис. 3 изображены рассчитан- ные отражательным методом сейсмограммы радиальных (R) и трансверсальных (T) компо- нент с источником, генерирующим qP-волну с сигналом Риккера частотой f=30 Гц. Расчет был выполнен для среды с трещиноватостью без поглощения (ΔN=0,6; ΔT=0,53) и для среды с по- глощающей трещиноватостью (ΔN=0,6+i0,054; ΔT=0,53+i0,004). Для данной модели среды была рассчитана и выведена в прямоугольной си- стеме координат (рис. 4, а — непрерывная линия) индикатриса групповых скоростей. Групповые скорости практически не зависят от наличия поглощения в трещиноватости. По Рис. 2. Схема наблюдений. Источник расположен в начале системы координат, приемники находятся на расстоянии 1 км от источника в вертикальной плоскости, перпендику- лярной плоскостям трещиноватости (R, T — направления регистрируемых колебаний вдоль радиуса окружности и вдоль касательной соответственно). Рис. 3. Сейсмограммы R- и T-компонент в зависимости от угла с вертикалью для трещиноватой среды без поглощения (а, б) и с поглощением (в, г). плоскости которой перпендикулярны оси OX. Свойства вмещающей среды определяются скоростью продольных волн VP=2800 м/с, по- РАСЧЕТ ВОЛНОВЫХ ПОЛЕЙ ДЛЯ АНИЗОТРОПНЫХ СРЕД С ПОГЛОЩЕНИЕМ НА ОСНОВЕ... Геофизический журнал № 4, Т. 33, 2011 125 Аки К., Ричардс П. Количественная сейсмология. Том 1. — Москва: Мир, 1983. — 519 с. Каппелини В., Константинидис А. Дж., Эмилиани П. Цифровые фильтры и их применение. — Мо- сква: Энергоатомиздат, 1983. — 360 с. Молотков Л. А. Матричный метод в теории распро- странения волн в слоистых упругих и жидких средах. — Ленинград: Наука, 1984. — 201 с. Роганов В. Ю. Вычисление волновых полей для ани- зотропных сред с использованием метода Ха- скелла—Томсона // Геофиз. журн. — 2009. — 31, № 3. — С. 63—73. Рис. 4. Значения групповых скоростей V (а — кружочки) и добротностей Q (б — кружочки), вычисленные по синтези- рованному волновому полю, и теоретические кривые (непрерывные линии). временам зарегистрированных сигналов также были найдены групповые скорости и нанесе- ны на построенный график кружочками (см. рис. 4, а). Из рисунка видно, что кинематиче- ские свойства сформированного волнового поля хорошо соответствуют теоретическим. Были рассчитаны теоретические значения добротностей Q в зависимости от группового угла α и представлены непрерывной кривой на рис. 4, б. Добротности также были вычис- лены на основе сравнения амплитуд и спек- тров зарегистрированных сигналов, распро- страняющихся в средах с поглощением и без поглощения. Результаты нанесены на график кружочками (см. рис. 4, б). Из сопоставления полученных значений Q с теоретическими значениями видно, что отра- жательный метод хорошо моделирует эффект анизотропного поглощения в трещиноватых средах. Выводы. 1. Разработан эффективный и устойчивый способ формирования 3D—3C волновых полей для многослойной анизотроп- ной горизонтально-слоистой среды на основе отражательного рекурсивного метода. Метод позволяет учитывать различные типы источ- ников, свободную поверхность, интерферен- ционные явления, связанные с тонкослои- стостью, наличие разно наклоненных систем трещин с поглощением. 2. Для градиентной среды выведено матрич- ное уравнение Риккати, коэффициенты кото- рого определяются по методу возмущений. 3. Результаты продемонстрированы на при- мере однородной среды с наличием трещино- ватости с поглощением. Список литературы Abo-Zena A. Dispersion function computations for unlimited frequency values // Geophys. J. Roy As- tronom. Soc. — 1979. — 58, № 1. — P. 91—105. Abramovici F. Diagnostic diagrams and transfer func- tions for oceanic wave guides // Bull. Seism. Soc. Amer. — 1968. — 58. — P. 427—456. Chichinina T., Obolentseva I., Gik L., Bobrov B., Ronquil- lo-Jarillo G. Attenuation anisotropy in the linear-slip model: Interpretation of physical modeling data// Geophysics. — 2009. — 74, № 5. — P. WB165— WB176. Dunkin I. W. Computation of modal solutions in layered Ю. В. РОГАНОВ, В.Ю. РОГАНОВ 126 Геофизический журнал № 4, Т. 33, 2011 elastic media at high frequencies // Bull. Seismol. Soc. Amer. — 1965. — 55, № 2. — P. 335—358. Fuchs K., Muller G. Computation of synthetic seismo- grams with the reflectivity method and comparison with observations // Geophys. J. Roy. Astronom. Soc. — 1971. — 23. — P. 417—433. Kennett B. L. N. Reflections, rays and reverberations // Bull. Seism. Soc. Amer. — 1974. — 64. — P. 1685— 1696. Kennett B. L. N. The effect of attenuation of seismo- grams // Bull. Seism. Soc. Amer. — 1975. — 65. — P. 1643—1651. Kennett B. L. N. Theoretical reflection seismograms for an elastic medium // Geophys. Prosp. — 1979. — 27, № 2. — P. 301—321. Kennett B. L. N., Kerry N. J. Seismic waves in stratified half-space // Geophys. J. Roy. Astronom. Soc. — 1979. — 57, № 3. — P. 557—583. Kennett B. L. N. Seismic wave propagation in stratified media. — New York: Cambridge University Press, 1983. — 200 p. Kind R. Computation of reflection coefficients for lay- ered media // J. Geophys. — 1976. — 42. — P. 191— 200. Fryer G. J. A slowness approach to the reflectivity me- thod of seismogram synthesis// Geophys. J. Roy. Astronom. Soc. — 1981. — 63. — P. 747—758. Fryer G. J., Frazer L. N. Seismic waves in stratified aniso- tropic media // Geophys. J. Roy. Astronom. Soc. — 1984. — 78. — P. 691—710. Mallick S., Frazer L. N. Practical aspects of reflectivity modeling // Geophysics. — 1987. — 52, № 10. — P. 1355—1364. Mallick S., Frazer L. N. Rapid computation of multioffset vertical seismic profile synthetic seismograms for layered media // Geophysics. — 1988. — 53, № 4. — P. 479—491. Norris A. N., Shuvalov A. L. Wave impedance matrices for cylindrically anisotropic radially inhomogeneous elastic solids // Quart. J. Mech. Appl. Math. — 2010. — 63, № 3. — P. 1—39. Sipkin S. A., Orcutt J. A., Jordan T. H. An examination of ScS travel times with causal Q reflectivity algorithm for SH polarized waves // Trans. Am. Geoph. Union — 1978. — 59. — P. 324. Stephen R. A. Synthetic seismograms for the case of the receiver within the reflectivity zone // Geophys. J. Roy. Astronom. Soc. — 1977. — 51. — P. 169—181. O’Neil M. E., Hill D. P. Causal absorption: its effect on synthetic seismograms computed by the reflecti- vity method // Bull. Seis. Soc. Am. — 1979. — 69. — P. 17—26. Ursin B. Review of elastic and electromagnetic wave propagation in horizontally layered media // Geo- physics. — 1983. — 48, № 8. — P. 1063—1081.