Акустический журнал, 2020, T. 66, № 2, стр. 198-205

ВВЕДЕНИЕ

Исследованный ранее авторами метод микросейсмического зондирования (ММЗ) [1–4] опирается на представление об определяющем вкладе в энергию микросейсмического поля фундаментальной моды волны Рэлея. Такое представление позволяет без дополнительных условий применять ММЗ для диапазона глубин земной коры с низким градиентом скоростей по глубине, то есть там, где амплитуды высших мод существенно ниже амплитуды фундаментальной моды. Однако, даже в таких благоприятных условиях применимость ММЗ зачастую вызывает вопросы ввиду важного обстоятельства. В ряде экспериментальных исследований показано, что в определенных условиях микросейсмическое поле может содержать заметную долю волн объемного типа. В частности, Л.П. Винник с соавторами показал, что в тихой внутриконтинентальной области доля объемных волн может составлять до 60% энергии микросейсмического поля в области частот около 1 Гц [5, 6]. Анализ функций взаимной корреляции записей микросейсмического фона в области 1 Гц, проведенный на территории штата Юта (Uinta Basin) США [7], свидетельствует о присутствии объемных волн в микросейсмическом поле от сильно удаленных слабых источников. Сравнительно недавние исследования также свидетельствуют о наличии в составе микросейсмического поля P-волн, источником которых являются удаленные штормовые процессы и ураганы [8–11]. Преобладание в микросейсмическом сигнале волн поверхностного или объемного типа, по всей видимости, определяется региональными особенностями самого явления [6].

В этой связи необходимо выяснить, как присутствие большой доли объемных волн в микросейсмическом поле искажает результат применения ММЗ. Правомерен вопрос, возможно ли вообще применять ММЗ (даже в благоприятных условиях слабого присутствия высших мод Рэлея), если в зондирующем микросейсмическом сигнале до 60% энергии формируется объемными волнами? Достаточно очевидно, что ввиду использования в ММЗ исключительно вертикальных компонент, а также с учетом уменьшения скоростей сейсмических волн при приближении к поверхности, в данном случае достаточно проанализировать осложнение микросейсмического поля за счет только объемных продольных сейсмических волн с падением, близким к вертикальному.

Важно отметить, что помимо поверхностных волн Рэлея и объемных волн на поверхности могут регистрироваться волны Лява. Однако, в ММЗ используются только вертикальные компоненты микросейсмического поля, поэтому ввиду своей природы волны Лява не участвуют ни в формировании сигнала, ни в его зашумлении.

Известно, что модель формирования сигнала в геофизических методах является основой построения процедуры инверсии – точного или приближенного решения обратной задачи. Ввиду сложности решения обратной задачи в общем виде, как правило, прибегают к упрощению модели формирования сигнала, снижая размерность задачи или используя то или иное допущение относительно свойств среды или свойств источника, либо используют более грубые физические закономерности в самой модели [12–14]. Проверка выполнения условий модели может требоваться в каждом конкретном случае.

В рамках решения поставленной задачи авторы исследовали закономерности формирования амплитудной реакции модельного микросейсмичeского поля на заглубленную неоднородность в присутствии поверхностных и объемных продольных волн на численной модели с использованием метода конечных разностей.

ОПИСАНИЕ МЕТОДА

Математическая модель, построенная из первых принципов на основе II-го закона Ньютона в дифференциальной форме для сплошной среды (1) и обобщенного закона Гука (2), использовалась для описания моделируемой системы в приближении небольших деформаций (3) [15, 16]:

где ${\mathbf{u}}({\mathbf{r}},t)$ – векторное поле смещений при упругой деформации среды, ${{\sigma }_{{ij}}}$ – тензор упругих напряжений материала среды, ${\hat {u}}$ – тензор деформации, $\rho ({\mathbf{r}})$ – плотность среды, $\lambda ({\mathbf{r}})$ и $\mu ({\mathbf{r}})$ – 1-й и 2-й упругие модули Ламэ (функции координат), ${\mathbf{f}}({\mathbf{r}},t)$ – плотность внешней силы, ${\hat {I}}$ – единичный тензор, r – радиус-вектор точки полупространства, t – время. Далее без ограничения общности будем рассматривать задачу в двумерной постановке, поскольку для описания исследуемых в рамках данной работы явлений не требуется большей размерности задачи. В расчетной области будем использовать прямоугольную Декартову систему координат xy, где y – вертикальная ось, направленная к центру Земли, граница y = 0 – дневная поверхность. Система уравнений (1)–(3) замыкается нулевыми начальными условиями ${{u}_{i}}({\mathbf{r}},t \leqslant 0) = 0,$ $\frac{{\partial {{u}_{i}}}}{{\partial t}}({\mathbf{r}},t \leqslant 0) = 0,$ а также следующими граничными условиями:

1) на свободной поверхности y = 0 условие отсутствия нормальных напряжений ${{\sigma }_{{xy}}} = {{\sigma }_{{yy}}} = 0$ приводит к системе уравнений

где $\nu (x,y) = \frac{\lambda }{{2(\lambda + \mu )}}$ – коэффициент Пуассона;

2) вертикальная граница расчетной области $x = 0$ – источник фундаментальной моды волны Рэлея заданной амплитуды A (по вертикали) на поверхности и циклической частоты ω_n [17]:

где A – амплитуда вертикальных колебаний на поверхности, φ – начальная фаза, $W\left( {t,\theta ~} \right)$ – огибающая с параметром нарастания θ, имеющая плавную форму, например, в виде функции Ханнинга

Вспомогательные константы q и s, зависящие от волновых чисел продольных, поперечных и рэлеевских волн k_P, k_S, k_R, во вмещающей среде

где значение волнового числа фундаментальной моды волны Рэлея k_R во вмещающей среде явно определяется с использованием точной формулы Малишевского [18]

3) нижняя горизонтальная граница y = Y – источник восходящей плоской P-волны единичной амплитуды с той же сеткой частот ω_n:

4) граничные условия (ГУ) на удаленной вертикальной границе в случае расчетов с волной Рэлея ${{{\Gamma }}_{\infty }}$ – отражающая граница:

5) периодическое граничное условие вдоль оси x в расчетах с объемной P-волной, что позволяет избежать возникновения краевых артефактов и разрыва фронта плоской волны:

Таким образом, получаем две задачи, отличающиеся набором граничных условий – первая система уравнений (1)–(3) с ГУ (4), (5), (8) для случая облучения волной Рэлея (см. рис. 1) и вторая система (1)–(3) с ГУ (4), (7), (9) для случая облучения исследуемой области восходящей продольной волной. Для численного интегрирования указанных уравнений используем в прямоугольной расчетной области xy равномерную квадратную сетку $0 \leqslant i < {{N}_{X}},$ $0 \leqslant j < {{N}_{Y}}$ с шагом, одинаковым по обеим осям. Для сеточного аналога векторного поля ${{u}_{k}}\left( {x,y,t} \right)$ смещений точек среды на введенной сетке используем обозначение $U_{k}^{{ijt}},$ где $k = x,y$ – компоненты вектора смещений.

Правую часть уравнения (1) без внешней силы f с учетом (2) и (3) можно записать в виде

где $i,~j,~k \in \overline {x,~y} $, ${{\delta }_{{ij}}}$ – дельта-символ Кронекера, и подразумевается соглашение Эйнштейна о суммировании по повторяющимся индексам.

Сеточная запись выражения (10) с использованием схемы второго порядка аппроксимации для всех производных по пространственным переменным (11), (12)

где в качестве w может быть любая функция от x, y, имеет следующий вид где производные сеточных функций ${{\lambda }^{{ij}}}$ и ${{\mu }^{{ij}}}$ по соответствующим переменным обозначены через $\lambda _{x}^{{ij}} = \frac{{\partial {{\lambda }^{{ij}}}}}{{\partial x}},$ $\mu _{x}^{{ij}} = \frac{{\partial {{\mu }^{{ij}}}}}{{\partial x}}$ и т.д.

Симметричная схема для второй производной по времени также имеет второй порядок аппроксимации

где для текущего временного шага индекс t опущен $U_{k}^{{ij}} = U_{k}^{{ijt}},$ а для прошлого и будущего временных шагов использованы диакритические знаки $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{U} _{k}^{{ij}} = U_{k}^{{ij(t - 1)}}$ и $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{U} _{k}^{{ij}} = U_{k}^{{ij(t + 1)}}$ соответственно, по аналогии с принятыми обозначениями в [19]. С учетом (13)–(16) и (17) для интегрирования основных уравнений (1)–(3) получаем явную схему во внутренних (не граничных) точках области:

Для интегрирования граничного условия на свободной поверхности (4) с использованием несимметричной сеточной разности второго порядка аппроксимации

можно построить неявную схему или явную в виде где $0 < i < {{N}_{X}} - 1$ и введен параметр $\vartheta = \frac{\nu }{{1 - \nu }},$ зависящий от коэффициента Пуассона. Система линейных алгебраических уравнений (СЛАУ) (19) или (20) относительно компонент вектора смещений в точках поверхности $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{U} _{k}^{{i0}},$ $0 < i < {{N}_{X}} - 1,$ $k = x,y$ дополняется тождествами, следующими из краевых условий типа Дирихле на границах $i = 0$ и $i = {{N}_{X}} - 1$ в случае облучения рэлеевской волной, либо периодическим граничным условием в случае облучения P-волной. Правые части уравнений в СЛАУ (19) и (20) содержат значения компонент вектора смещения $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{U} _{k}^{{i1}}$ и $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{U} _{k}^{{i2}}$ в двойном подповерхностном слое j = 1, 2 на будущем временном шаге t + 1, вычисляемые явным образом на каждом текущем временном шаге t из системы (18). Таким образом, обе СЛАУ (19) и (20) с учетом граничных условий и значений во внутренних точках $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{U} _{k}^{{i1}}$ и $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{U} _{k}^{{i2}},$ $0 < i < {{N}_{X}} - 1$, $k = x,y$ представляют собой полные системы.

Шаг интегрирования по времени τ выбирается согласно теореме Котельникова в зависимости от максимальной частоты ω_n колебаний в системе $\tau \leqslant \frac{1}{\zeta }\frac{\pi }{{{{\omega }_{n}}}},$ а также из условия сходимости схемы согласно обобщенному условию Куранта [19], определяемому выражением $\tau = \frac{1}{\chi }\frac{h}{{{{V}_{{P\max }}}\sqrt 2 }},$ где ξ, $\chi \geqslant 1$ – коэффициенты запаса прочности, ${{V}_{{P\max }}}$ – максимальное значение скорости распространения P-волн в модельной среде. Реализацию описанного подхода для построения модели в трехмерном случае, а также описание численного алгоритма и его верификацию можно найти в [1] или более подробно в [20].

Физические параметры среды и неоднородности выбраны в соответствии с вариантами, рассмотренными в предыдущих исследованиях [1, 2, 21, 22] для возможности сопоставления получаемых результатов моделирования:

1) параметры вмещающей среды: v_P = 3193.74 м/с, v_S = 1843.91 м/с, v_R = 1695.29 м/с, упругие модули Ламэ λ = μ = 6.8 ГПа;

2) низкоскоростная неоднородность: v_P = = 2554.99 м/с, v_S = 1475.13 м/с, λ = μ = 4.352 ГПа;

3) высокоскоростная неоднородность: v_P = = 3992.17 м/с, v_S = 2304.89 м/с, λ = μ = 10.625 ГПа.

Плотности материала вмещающей среды и неоднородности одинаковы и составляют ρ = 2000 кг/м³, коэффициент Пуассона также сохраняет значение по всему объему ν = 0.25. Горизонтальный размер рассмотренных неоднородностей L = 2, 10, 20 км, вертикальный d = 1 км, глубина залегания кровли h₁ = 1.5 км, подошвы h₂ = 2.5 км (см. схему на рис. 1). Размер расчетной области 80 × 36 км при шаге сетки h = 100 м, всего 288 000 ячеек квадратной сетки. Шаг интегрирования по времени составлял τ = 7.83 мс, обеспечивающий сходимость интегрирования. Продолжительность каждого расчета составляла 6028 временных шагов, что соответствует 47.2 с модельного времени.

Набор частот f_n гармонических источников P- и R-волн выбирался в полосе 0.125–1.695 Гц таким образом, чтобы обеспечить равномерное покрытие диапазона длин волн фундаментальной моды волны Рэлея λ_R(f_n) = 1000, …, 13600 м с шагом 200 м (всего 64 частоты). С использованием параллельных вычислений для каждого численного эксперимента в рамках одного расчета на 64 ядрах выполнялось облучение 64 волнами разной частоты. При этом на каждом расчетном ядре облучение происходило независимо, после чего данные собирались первым процессом и подготавливались к анализу.

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

С использованием описанных моделей была проведена серия из 12 численных экспериментов – для двух противоположных скоростных контрастов неоднородности (0.8 – низкоскоростная по отношению к вмещающему полупространству, 1.25 – высокоскоростная), для трех различных латеральных размеров (2 км – меньше длины волны Рэлея ${\lambda }_{R}^{*}$ на частоте максимальной амплитудной реакции f *, 10 км – сравнима с длиной волны ${\lambda }_{R}^{*}$ и 20 км – больше нее) и для двух типов источника (объемные продольные волны с вертикальным направлением прихода или волны Рэлея, облучающие сбоку). Важно отметить, что вертикальный размер неоднородности во всех численных экспериментах составлял 1 км, что составляет 1/8–1/4 реагирующей длины волны Рэлея.

Синтетические изображения модельных неоднородностей, полученные в соответствии с процедурой ММЗ, примененной к синтетическому полю на поверхности, которое сформировано в результате рассеяния P-волн, R-волн и равнозначной смеси из P- и R-волн, представлены на рис. 2 и 3. Изображения строились путем вычисления квадрата амплитуды вертикальной компоненты смещения в каждой точке поверхности как функции частоты и расстояния вдоль профиля $V\left( {x,f} \right)$. Облучение производилось волнами с единичной амплитудой вертикальных колебаний. Согласно процедуре ММЗ из пространства (x, f) осуществлялся переход в пространство (x, λR(f)) или (x, h), где h – глубина, определяемая соотношением

k_G – коэффициент глубинной привязки [1].

Рассмотрим вначале случай с низкоскоростной неоднородностью, малой по сравнению с длиной падающей волны (рис.2). Видно, что над основной аномалией, формирующей изображение исходного объекта, присутствуют значительные по величине аномалии того же знака, что и в основной (рис. 2а). При облучении волной Рэлея основная аномалия проявляется в тех же частотах, что и в продольной волне, но над основным изображением формируется аномалия обратного знака по отношению к основной (рис. 2г). Когда модельное зондирующее поле представлено обоими типами волн в равных долях (рис. 2ж), основная аномалия имеет размер и форму, близкие к случаям отдельных облучений (рис. 2а и 2г), однако артефакты в зоне коротких волн над основным изображением взаимно компенсируют друг друга. В этом случае можно даже говорить о своеобразном улучшении изображения при использовании ММЗ, если около 50% энергии микросейсмического сигнала составляют P-волны.

При облучении P-волной низкоскоростной неоднородности с латеральным размером, сравнимым с длиной падающей волны (рис. 2б, 2д, 2з), основная аномалия имеет более искаженную форму границы (рис. 2б), чем при облучении R-волной (рис. 2д). Как и в случае малой неоднородности, формируется ложная положительная аномалия в коротких волнах. При облучении волнами обоих типов (рис. 2з) аномалия взаимно компенсируется. Присутствие 50% P-волн в облучающем поле (рис. 2з) искажает форму изображения по сравнению с облучением чистой волной Рэлея (рис. 2д). Незначительно увеличивается глубина изображения неоднородности. Аналогичные эффекты наблюдаются в случае протяженной по горизонтали неоднородности (рис. 2в, 2е, 2и).

Результаты моделирования с высокоскоростной неоднородностью представлены на рис. 3. Как видно, высокоскоростное включение проявляется падением интенсивности в определенной области спектра, связанной с глубиной и вертикальным размером включения. Изображение малой неоднородности по ММЗ в поле объемных волн (рис. 3а) близко к изображению в поле поверхностных волн (рис. 3г). Различие состоит в присутствии малой интенсивной аномалии в коротких волнах. При облучении P-волной основная аномалия более четкая и имеет несколько меньший горизонтальный размер, что более контрастно выделяет объект (рис. 3а). Результат при смешанном облучении (50% P-волны, 50% R-волны) малой неоднородности (рис. 3ж) в сравнении с вариантом рис. 3г имеет более интенсивный “обратный эффект” над неоднородностью, а также более четкое основное изображение неоднородности. Оба варианта с крупным высокоскоростным включением (рис. 3б, 3д, 3з и рис. 3в, 3е, 3и) имеют сходные черты. В присутствии P-волн усиливается интенсивность “обратного эффекта” в области коротких длин волн над неоднородностью. При этом лучше выделяется положение субгоризонтальных границ неоднородности, верхняя граница приобретает более ровную форму, нижняя граница становиться более контрастной, чем при облучении только волнами рэлеевского типа.

В обоих рассмотренных вариантах основная аномалия, интерпретируемая в ММЗ как искомый объект, имеет близкие к заданной скоростной неоднородности форму, “знак” и размеры. Она проявляется в одном и том же частотном диапазоне как в поле только вертикально падающих P-волн, так и при облучении только R-волнами. Этот результат сохраняется независимо от горизонтального размера и характера неоднородности. Видно, что любое содержание объемных P-волн в модельном микросейсмическом сигнале в основном сохраняет главную аномалию, но при этом могут возникать артефакты, не связанные со строением среды.

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

ЗАКЛЮЧЕНИЕ

Результаты моделирования на основе первых принципов для взаимодействия упругих P- и R-волн в линейном изотропном полупространстве с заглубленной скоростной неоднородностью позволяют в дальнейшем полагать, что в случае, когда вертикальная компонента микросейсмического поля помимо волн типа Рэлея содержит объемные продольные волны, использование метода микросейсмического зондирования для оценки строения геологической среды допустимо. Даже присутствие объемных P-волн вплоть до 100% по энергии не ограничивает возможности применения ММЗ без каких-либо процедурных изменений для оценки структуры среды. При этом ожидаемый результат может осложняться более интенсивными искусственными аномалиями в высокочастотной области, чем в случае, когда микросейсмическое поле состоит только из волн Рэлея. Оценка по ММЗ скорости в обнаруженном включении также будет искажена по сравнению с чисто рэлеевским полем. Это следует принимать во внимание при геологической интерпретации результатов ММЗ.

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

Работа выполнена с использованием оборудования Центра коллективного пользования сверхвысокопроизводительными вычислительными ресурсами МГУ имени М.В. Ломоносова [23, 24].