Математическое моделирование эффективных упругих параметров Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

УДК 519.677:550.834.01 Б01: 10.14529/ттр180201

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ЭФФЕКТИВНЫХ УПРУГИХ ПАРАМЕТРОВ

П.Н. Александров1, В. Н. Кризскии2

1Институт физики Земли им. О.Ю. Шмидта РАН, г. Москва, Российская Федерация 2Стерлитамакский филиал Башкирского государственного университета, г. Стерлитамак, Российская Федерация

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

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

Введение

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

Первый подход можно подразделить на две части - первая связана с аппроксимацией микронеоднородной упругой среды, заданной моделью вида материальных уравнений (законом Рука), которая требует определения окрестности компакта, т.е. вида материальных уравнений. Этот подход развивается И.О. Баюк и изложен в ее докторской диссертации [1]. Другая часть связана с осреднением закона Рука микронеоднородной упругой среды с учетом длины сейсмической волны [2-4]. В последнем случае, для наиболее простой одномерной тонкослоистой модели среды, пачки переслаивающихся пропластков, получен осредненный закон Рука, который провидит к

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

В общем случае, связь напряжений и деформаций в законе Гука, которая в векторном представлении будет иметь вид Р = Не, где Р - 9-ти компонентный вектор напряжений, е - 9-ти компонентный вектор деформаций, является линейной и описывается матрицей упругих параметров Н размерности 9x9 элементов, которая может зависеть от временной частоты в частотной области и, таким образом, описывать дисперсию упругих параметров.

1. Решение прямой задачи сейсморазведки для слоисто-анизотропных сред

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

Рассмотрим модель горизонтально-слоистой анизотропной геологической среды в декартовой системе координат с осью г, направленной вглубь Земли. Вводя вектор

действием вертикальных сил, Я - вектор смещений, и применяя преобразование Фурье по горизонтальным координатам х, у и по времени систему уравнений Ламе для этой модели можно свести к виду

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

Продолжая поле X через горизонтально слоистую среду от дневной поверхности

j-ro слоя, п - номер последнего слоя бесконечной толщины, X0 - вектор-столбец, заданный на дневной поверхности, zn - глубина залегания последней границы, Aj -передаточная матрица j-ro слоя. Поскольку плоскость z = z0 совпадает с поверхностью Земля/Воздух, то выполняется условие X*(z0) = X0.

В нижнем слое бесконечной толщины^ основываясь на знаке действительной ча~ сти собственных значений матрицы An, го представления Xn(z) = e (z Zn)Xn(zn) выделим решение X— возрастающее при z ^ и решение X+, убывающее при z ^ Xn(z) = X+ + X-. Удовлетворяя условию на бесконечности, необходимо

положить X+ = CSC~lXn(z) = 0 везде, в том числе и при z ^ zn. Здесь C - матрица, составленная из собственных

An SS

компонентный вектор напряжении, появляющийся под

до подошвы последнего слоя, получим Xn(zn)

CSC

Z = — d11 1d

42

—d21 ld-

■22-

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

Отсюда вытекает связь между компонента ми поля X0: C SC -1Xn(z) =

и eAjМ X0 = DX0 = 0, D =( d11 d12 j=0 I V d21 d22

На дневной поверхности зададим сторонние напряжения Pz = Pf (например, в виде плиты вибратора конечных (реальных) размеров), тогд& решение поставленной задачи будет иметь вид S = Z-1Pszt. При этом никаких ограничений на структуру матрицы упругих параметров не накладывается, и, кроме этого, она может быть частотно зависимая.

Используя обратное преобразование Фурье по пространственным и временной частотам, можно перейти в пространственно-временную область.

Вычислительные эксперименты проведены для модели трехслойной среды (рис. 1), у которой второй слой является макроанизотропньтм и эффективные параметры которого являются результатом осреднения тонкослоистой среды с параметрами прослоек указанных в таблице. В таблице использованы следующие обозначения: h - толщина пропластка, Vp - скорость продольных волн в слое, Vs - скорость поперечных волн в слое, р - плотность горной породы в слое.

Рис. 1. Модель трехслойной среды с квазианизотропным вторым слоем. Параметры квазианизотропии вычислены в результате осреднения тонкослоистой среды, с наклоном слоистости под углом 45° к ос и у

Таблица

Параметры тонкослоистой пачки

h, м Vp, м/с Vs, м/с р, кг/м3

0,01 565,6854 332,7561 2000

0,02 126,4911 74,4065 2000

0,03 565,6854 332,7561 2000

0,01 126,4911 74,4065 2000

0,02 565,6854 332,7561 2000

Для каждой прослойки второго слоя вычислялись параметры Ламе по формулам: /р2

X = p(Vp2 — 2Vs2), л = pVs2- После чего проводилось осреднение и осуществлялся

1 0 0

поворот с направляющими косинусами | 0 0, 7071 —0, 7071

0 0, 7071 0, 7071

В случае воздействия на среду импульсом сторонних вертикальных напряжений, заданных по площадке 1x1 м в начале координат с формой, представленной на рис. 2, вычислительным экспериментом получено сейсмическое поле, фрагмент вертикальной компоненты которого изображен на рис. 3. Здесь: момент времени Ь = 0, 31 с, площадка измерений - квадрат на плоскости г = 0 размерами 120x120 м с центром в начале координат.

Форма импульса

■0 1 -0.08 -О 0G -0 04 -0.02 0 0.02 0 04 0.06 0 0В 0 1

Время (сек.)

Рис. 2. Форма импульса для сторонних вертикальных напряжений, заданных по площадке 1x1 м в начале координат

Sz (t = 0.31 сек)

■60 -40 -20 0 20 40 60

Y, м

Рис. 3. Фрагмент вертикальной компоненты сейсмического поля, рассчитанного в момент времени t = 0, 31 с по площади 120x120 м

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

2. Соответствие эффективных сейсмических скоростей их средневзвешнным величинам

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

В однородном анизотропном пространстве уравнения теории упругости можно записать в виде

д2 1 <й2Ч = - м Ч,

дЬ2 р

где м = В * Н * Е, В =

гкх гку гкх 0 0 0 0 0 0

0 0 0 %кх гку гкх 0 0 0

0 0 0 0 0 0 гкх гку гк

Н

г,з = 1, 9 Е

гкх [1] %ку [1] гкх [1]

Для случая плоской волны, распространяющейся в направлении оси г, имеем

%кх = гку = 0, отсюда

8 = - к 21

о, 2 8 — кх

дЬ2 р

Ух

< ± >-1 М

0

< 1 >-1

А+2м

0

М

00

122,878 0

0 122,878 0 0 0 208,8932

Б,

< 0

>

-1

Для случая плоской волны, распространяющейся в направлении оси у, имеем %кх = гкх = 0, отсюда

д2 1 —Ч = - к о1

дЬ2 Ч ку р

0

0 < Щ+М > + < А+М

< ц>

0 0

>2 <

А+2М

>

1

0

0 0

< 1 >"1

Ч,

/ 275,0695 0 0

Уу = I 0 449,5187 0 \ 0 0 122,878

Для случая плоской волны, распространяющейся в направлении оси х, имеем

%ку = %кх = 0, отсюда

Ч = _ к о 1

гл,2 Ч = кх

дЬ2 р

< 4м(а+М) > + < А >2 < 1 >-1

< А+2м > + < А+2м > < А+2м >

0 0

0

< Ц>

0 0

0 < 1 >-1

Ч,

х

1

/449,5187 0 0

К = I 0 275,0695 0 \ 0 0 122,878

Анализируя коэффициент волнового уравнения, который характеризует матери-эльные свойства вещества (горной породы), получены скорости, количество которых равно 9, из которых 5 равны друг другу. В результате, в анизотропной среде распространяются волны с разными скоростями в разных направлениях. Для рассматриваемой модели независимых скоростей четыре: Ур1 = 449,5187 м/с, = 275,0695 м/с, Ур2 = 208, 8932 м/с, У,2 = 122, 878 м/с.

С другой стороны, можно вычислить скорости, основываясь не на динамическом подходе (через уравнения Ламе), а на кинематическом подходе (через средние скорости и средние интервальные скорости). Для средних скоростей, вычисленных на основе уравнения эйконала, получим: Е ЫУ;

< Vp>i

< У>

5

Е hi

i=1

Е hi Vi

i=1 5

hi

i=1

419,2873 м/с

< Vp>m

1

5 h-V hhx

■Л vi i=1 vp

hi

i=1

246,6396 м/с ,<V4 >m

i

5 hi i=1 vi

5

hi

i=1

5

hi

i=1 5 и hi vi i=1 vp

5

hi

i=1 _

5 и ' hi

i=1 vi

262, 2096 м/с,

154, 2409 м/с.

Эти величины скоростей не совпадают со скоростями, полученными по осреднению закона Гука.

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

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

3. Уравнения эйконала в произвольно-анизотропной неоднородной среде

В произвольно-анизотропной и неоднородной среде уравнения эйконала являются собственными значениями матрицы размерности 3x3 [5]:

A

т

gradt [0] [0] [0] gradt [0] I H [0] [0] gradt

gradt [0] [0] [0] gradt [0] [0] [0] gradt

1 0 0 010 001

Следовательно, в произвольно-анизотропной и неоднородной среде имеется три уравнения эйконала вида

У gradí = 1.

P

Матрицу упругих параметров можно представить в виде составной матрицы

к11 к12 к

Н= к21 к22 к

к31 к32 к

13 23 33

где подматрицы к^ являются матрицами размерности 3x3. При этом, исходя из структуры матрицы Л, в построении уравнений эйконала будут участвовать только диагональные подматрицы к11, к22, к33. Остальные подматрицы матрицы упругих параметров не участвуют в формировании скорости распространения упругого поля.

Рассмотрим решение уравнения эйконала. Пусть матрица связана с одной из под-к11 к22 к33. Решение уравнения эйконала вида Уgradí = 1, где

матрица У размерности 3x3 представима в виде У = V [А]V-1. Здесь матрица V есть матрица, составленная из собственных У

А11 0 0

[А] = 0 А22 0

0 0 А33

есть матрица собственных У

Пусть в однородной среде имеет место соотношение gradí = vgradí, тогда функция Ь будет удовлетворять уравнению эйконала вида

Или

Откуда

ГТ1 -""/ТТ 1 л- л-

gradí Уgradí = gradí V- Уvgradt = gradí [АJgradí = 1.

А"< |)2 + А22 (|)2 + А22( I)2 = 1

/

Ь = ^ х2 + ~у2 + г 2.

А11 А22 А33

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

Заключение

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

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

Литература

1. Баюк, И.О. Междисциплинарный подход к прогнозированию макроскопических и фильтрационно-емкостных свойств коллекторов углеводородов: дис. ... доктора физ.-мат. наук / И.О. Баюк. - М., 2013. - 228 с.

2. Уайт, Дж.Э. Возбуждение и распространение сейсмических волн / Дж.Э. Уайт. - М.: Недра, 1986.

3. Backus, G.E. Long-Wave Anisotropy Produced by Horizontal Layering / G.E. Backus // Journal of Geophysical Research. - 1962. - V. 67, № 11. - P. 4427-4440.

4. Рытов, C.M. Акустические свойства мелкослоистой среды / С.М. Рытов // Акустический журнал. - 1956. - Т. 2, № 1. - С. 71-83.

5. Александров, П.Н. Вывод уравнения эйконала для анизотропных неоднородных сред / П.Н. Александров // Юбилейная X ежегодная международная конференция «Гальпе-ринские чтения - 2010>. Инновационные технологии и фундаментальные исследования в наземно-скважинной сейсморазведке 2D, 3D, ВСП и сейсмологии, посвященная 90-летию Е.И. Гальперина. - М., 2010. - С. 53-59.

Павел Николаевич Александров, доктор физико-математических наук, главный научный сотрудник, ЦГЭМИ, Институт физики Земли им. О.Ю. Шмидта FAH (г. Москва, Госсийская Федерация), Alexandr@igemi.troitsk.ru.

Владимир Николаевич Кризский, доктор физико-математических наук, профессор, кафедра «Математическое моделирование>, Стерлитамакский филиал Башкирского государственного университета (г. Стерлитамак, Госсийская Федерация), Krizsky@rambler.ru.

Поступила в редакцию 24 а,прем,я 2018 г.

MSC 65Z05 DOI: 10.14529/mmpl80201

MATHEMATICAL MODELLING OF EFFECTIVE ELASTIC

PARAMETERS

P.N. Aleksandrov1, V.N. Krizsky2

Schmidt Institute of Fhysics of the Earth of the RAS, Moscow, Russian Federation,

2

E-mail: Alexandr@igemi.troitsk.ru, Krizsky@rambler.ru

The article is devoted to the study of the laws of elastic field propagation in inhomogeneous anisotropic media. At the same time, anisotropy is introduced as effective (averaged) parameters of a thin-layered medium, which determines the macroanisotropic elastic parameters of the rock. It is shown that the effective elastic parameters obtained from the theory of elasticity (Lame equations) do not coincide with the effective parameters

obtained using the kinematic approach. On the basis of reduction of equations of the theory of elasticity to the systems of the ordinary differential equations of the first order the solution of a direct problem of seismic exploration (as a boundary value problem) for horizontally layered and anisotropic model of the geological environment is received. The given result of calculation of the seismic field registered on the daily surface in the case of an anisotropic object leads to the complex picture of the wave field. This means that it is necessary to improve the methods of seismic exploration in the course of the studying of anisotropic properties of the geological environment.

Keywords: seismic anisotropy; effective parameters; systems of differential equations; eikonal equation.

References

1. Bayuk I.O. Mezhdisciplinarny podkhod k prognozirovaniyu makroskopicheskikh i filtraciormo emkostnykh svoystv kollektorov uglevodorodov [Interdisciplinary Approach to Prediction of Macroscopic and Filtration-Capacitive Properties of Hydrocarbon Reservoirs. Thesis of Doctor of Physical and Mathematical Sciences], Moscow, 2013. 228 p. (in Russian)

2. Wait J.E. Underground Sound - Application of Seismic Waves, Elsevier Science Publishers, Amsterdam, 1983.

3. Backus G.E. Long-Wave Anisotropy Produced by Horizontal Layering. Journal of Geophysical Research, 1962, vol. 67, no. 11, pp. 4427-4440. DOI: 10.1029/JZ067i011p04427

4. Rytov S.M. [Acoustic Properties of Fine-Layered Medium]. Akusticheskiy zhurnal, 1956, vol. 2, no. 1, pp. 71-83. (in Russian)

5. Aleksandrov P.N. [Derivation of the Eikonal Equation for Aanisotropic Inhomogeneous Media]. X Annual International Conference "Halperin Readings - 2010". Innovative Technologies and Fundamental Researches in 2Q, 3D VSP Well-Surface Exploration and Seismology Commemorating 90 Year Anniversary of Y.I. Galperin, Moscow, 2010, pp. 53-59. (in Russian)

Received April 24, 2018