Amplitude characteristics of Rayleigh-type waves in horizontally inhomogeneous layered media

Введение

Специфические свойства поверхностных акустических волн позволяют использовать их во многих практических приложения, например, в сейсморазведке [1-5] или дефектоскопии [6, 7]. В сейсморазведке они используются как средство исследования внутреннего строения внешних слоев Земли (определение границ слоев, дефектов и неоднородностей), поскольку скорость, затухание, амплитуды и структура рэлеевской волны связаны с упругими свойствами сред, вдоль границы которых она распространяется. Такие геофизические методы как MASW [1] и SPAC [2] используют информацию о скорости волны, и с их помощью восстанавливается модель исследуемой среды на основе анализа дисперсионной кривой. Другие подходы, например, метод микросейсмического зондирования (ММЗ) [3, 8] или метод спектрального отношения горизонтальной к вертикальной компоненте (H/V) [4, 9] основаны на измерении амплитуд и последующем анализе спектров измеренных сигналов. Оба подхода имеют свои сильные и слабые стороны [10] и возникает вопрос о возможности объединения этих методов для максимизации используемой экспериментальной информации и получения наиболее точного решения обратных задач. Для этого необходимо отчетливо понимать механизмы распространения поверхностной волны не только в вертикально, но и горизонтально неоднородных средах, причем, исследовать не только скоростные характеристики, но и амплитудные. Стоит отдельно отметить, что касательно этого вопроса были получены серьезные математические результаты [11]. Однако в силу сложности постановки и полученных аналитических решений их использование для прикладных задач существенно затрудняется. Настоящая работа нацелена на поиск наиболее простого и экономного подхода, но отражающего основные закономерности распространения ПАВ в горизонтально-неоднородных средах.

Исследованию свойств поверхностных волн в слоистых средах посвящен ряд работ, например, [12-18]. Однако преимущественно они посвящены анализу дисперсионных закономерностей, поэтому поведение амплитудных характеристик в различных неоднородных средах все еще нуждается в дальнейшем исследовании, поскольку они также несут информацию о строении среды [19], что актуально с прикладной точки зрения. Одной из таких характеристик является эллиптичность — отношение горизонтальной компоненты смещений к вертикальной, которая удобна для практических задач [4, 9], поскольку не требует учета особенностей источника, в отличие от отдельного исследования амплитуд компонент смещений, которое представляет собой существенно более сложную задачу [20]. На практике эта информация часто неизвестна, поэтому необходимо разрабатывать подходы, позволяющие описывать изменение амплитуды в отсутствие источника. Одной из таких моделей является модель плавного перехода между областями с различными значениями параметров среды [21]. Такая постановка позволяет определить амплитуду волны во второй среде по отношению к первой.

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

Постановка задачи

Рассматривается двумерный случай со следующей геометрией (рис. 1): свободная граница вертикально слоистой среды расположена при $z=0$, ось z направлена вниз. Слоистое упругое полупространство занимает область $z\ge 0$. Сверху при $-\infty <z<0$ находится вакуум. Вдоль свободной границы в положительном направлении оси x распространятся плоская поверхностная волна. Рассматриваемую геометрию разделим на три области: I - опорная среда, II - исследуемая среда, III - область перехода. Волна из области I переходит в область II, испытывая трансформацию в области III, в которой происходит плавное изменение параметров системы. В аналитических построениях волновое поле в областях I и II рассчитывалось независимо друг от друга. Связь между рассчитанными полями устанавливалась, исходя из предположения, что поток энергии через полубесконечное сечение сохраняется [21]. В численном моделировании были рассмотрены два случая: при длине области перехода равной одной длине волны и пятидесяти длинам волн.

Рис. 1. Постановка задачи.

Пав в плоско-слоистой среде

Выражения для акустического поля

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

Поле в каждом из однородных слоев описывается уравнением Ламэ, записанным для смещений в отсутствие источника [22-24]:

$\rho \frac{{\partial }^{2}u}{\partial t^2}=\mu \Delta u+(\lambda +\mu )\nabla \text{div}u$.                                      (1)

Смещение $u$ элемента упругой среды выражается через скалярный φ и векторный ψ потенциалы. При этом для двумерной задачи (в которой смещение вдоль координаты y отсутствует, а потенциалы от y не зависят) ψ = (0, ψ, 0), а смещения по осям x и z имеют вид:

 u =$\nabla$$\phi$ + rotψ,

 $u_x=\frac{\partial \phi }{\partial x}-\frac{\partial \psi }{\partial z}, u_z=\frac{\partial \phi }{\partial z}+\frac{\partial \psi }{\partial x}.$          (2)

Подставляя выражения (2) в (1) и применяя поочередно операции дивергенции и ротора, получим волновые уравнения для потенциалов

$\rho \frac{{\partial }^2\phi }{\partial t^2}-(\lambda +2\mu )\Delta \phi =\mathrm{0,}\rho \frac{{\partial }^2\psi }{\partial t^2}-\mu \Delta \psi =0$.

При гармонической зависимости от времени, они переходят в уравнения Гельмгольца:

$\begin{array}{l}\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial z^2}\right)\phi +k_l^2\phi =\mathrm{0,}\\ \left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial z^2}\right)\psi +k_t^2\psi =0.\end{array}$ (3)

В уравнениях (3) введены волновые числа $k_l=\omega /c_l$ - для продольных волн и $k_t=\omega /c_t$ - для поперечных. Решения данных уравнений отыскиваются в виде плоских волн:

$\phi =G\left(z\right)e^{-i\left(\omega t-kx\right)},\psi =F\left(z\right)e^{-i\left(\omega t-kx\right)}$,

где k = ω/с, с - фазовая скорость ПАВ.

Подставляя эти выражения в уравнения (3), получаем для функций $G\left(z\right)$ и $F\left(z\right)$ уравнения:

$\begin{array}{c}\frac{{\partial }^{2}F}{\partial z^2}-q^{2}F=\mathrm{0,}{q}^2=k^2-k_l^2,\frac{{\partial }^{2}G}{\partial z^2}-s^{2}G=\mathrm{0,} \\ s^2=k^2-k_t^2.\end{array}$

и их решения:

$F(z)=A^{(-)}{e}^{-qz}+A^{(+)}{e}^{qz},G(z)=B^{(-)}{e}^{-sz}+B^{(+)}{e}^{sz},$

где A^(-), B^(-), A⁽⁺⁾, B⁽⁺⁾ соответствуют амплитудам волн, распространяющихся в положительном и отрицательном направлении вдоль вертикальной оси. Отметим, что выражения для функций $G\left(z\right)$ и $F\left(z\right)$ являются общими решениями. Для записи решений в конкретном слое в дальнейшем будет использоваться соответствующий нижний индекс.

Так как при $z\to \infty$ поле должно затухать, то физический смысл имеют только решения, распространяющиеся в положительном направлении оси z. Соответственно, выражения для потенциалов в полупространстве (при $z>z_N$ ) будут выглядеть следующим образом:

$\phi =A_N{e}^{-q_{N}z}{e}^{-i\left(\omega t-kx\right)},\psi =B_N{e}^{-s_{N}z}{e}^{-i\left(\omega t-kx\right)}$.             (4)

Здесь величины $q_N=\sqrt{k^2-k_{lN}^2}$, $s_N=\sqrt{k^2-k_{tN}^2}$ определяют характерные масштабы затухания продольной и сдвиговой компонент поля вглубь полупространства. А для слоев решения представляются в следующем виде:

 $\begin{array}{l}{\phi }_i=A_i^{(-)}{e}^{-q_{i}z}{e}^{-i\left(\omega t-kx\right)}+A_i^{(+)}{e}^{q_{i}z}{e}^{-i\left(\omega t-kx\right)},\\ {\psi }_i=B_i^{(-)}{e}^{-s_{i}z}{e}^{-i\left(\omega t-kx\right)}+B_i^{(+)}{e}^{s_{i}z}{e}^{-i\left(\omega t-kx\right)},\end{array}$         (5)

где $q_i=\sqrt{k^2-k_{li}^2}$, $s_i=\sqrt{k^2-k_{ti}^2}$, i - номер слоя.

Граничные условия

Для определения характеристик ПАВ необходимо записать граничные условия, позволяющие найти скорость волны и связь между амплитудами потенциалов. Поскольку в системе из N слоев существует $4(N-1)+2$ волновых компонент, то и граничных условий должно быть $4(N-1)+2$. На свободной границе должно выполняться условие равенства нулю нормальных и касательных напряжений:

${\left.T_{zz}\right|}_{z=0}=\mathrm{0,}{\left.T_{xz}\right|}_{z=0}=0$.

На границе i-го и (i+1)-го слоя должно выполняться условие равенства нормальных и касательных смещений и напряжений:

$\begin{array}{c}{\left.u_z^i\right|}_{z=z_i}={\left.u_z^{i+1}\right|}_{z=z_i},{\left.u_x^i\right|}_{z=z_i}={\left.u_x^{i+1}\right|}_{z=z_i},\\ {\left.T_{zz}^i\right|}_{z=z_i}={\left.T_{zz}^{i+1}\right|}_{z=z_i},{\left.T_{xz}^i\right|}_{z=z_i}={\left.T_{xz}^{i+1}\right|}_{z=z_i},\end{array}$

где

$T_{zz}=\lambda \left(\frac{{\partial }^2\phi }{\partial x^2}+\frac{{\partial }^2\psi }{\partial z^2}\right)+2\mu \left(\frac{{\partial }^2\phi }{\partial z^2}+\frac{{\partial }^2\psi }{\partial x\partial z}\right),$ 

$T_{xz}=\mu \left(2\frac{{\partial }^2\phi }{\partial x\partial z}-\frac{{\partial }^2\psi }{\partial z^2}+\frac{{\partial }^2\psi }{\partial x^2}\right)$.

Характеристическая матрица

Подставляя решения в виде (4), (5) в граничные условия, получаем систему линейных уравнений относительно неизвестных амплитуд потенциалов, которую можно записать с использованием характеристической матрицы размерности $4(N-1)+2$. Равенство нулю определителя этой матрицы представляет собой дисперсионное уравнение, корнем которого является скорость ПАВ.

Условием существования нетривиального решения однородной системы линейных алгебраических уравнений является равенство нулю ее определителя. Так как при подстановке скорости ПАВ в матрицу, ее определитель обратится в ноль, то уравнения станут линейно зависимы и, следовательно, одно из них можно опустить. Затем, перенося один из столбцов в правую часть, получаем систему из $4(N-1)+1$ уравнений с отличной от нуля правой частью. Для поиска решения данной системы удобно применить метод Крамера [25]. Найденные значения являются амплитудными коэффициентами, связывающими амплитуды потенциалов с амплитудой потенциала, соответствующего перенесенному в правую часть столбцу. Эта связь позволяет записать выражения для потенциалов и, следовательно, для смещений через единственную нормировочную амплитуду.

Продемонстрируем нахождение амплитудных коэффициентов на примере двухслойной среды (т.е. слоя, лежащего на полупространстве). В этом случае характеристическая матрица и соответствующее дисперсионное уравнение выглядят следующим образом [22]:

 $d=\left(\begin{array}{cccccc}0& 0& k^2{\lambda }_1-{\alpha }_1{q}_1^2& -k^2{\lambda }_1+{\alpha }_1{q}_1^2& -2ik{\mu }_1{s}_1& 2ik{\mu }_1{s}_1\\ 0& 0& -2ik{\mu }_1{q}_1& 2ik{\mu }_1{q}_1& -{\mu }_1(k^2+s_1^2)& -{\mu }_1(k^2+s_1^2)\\ ik& s_2& -e^{-q_{1}h}ik& -e^{q_{1}h}ik& -e^{-s_{1}h}{s}_1& e^{s_{1}h}{s}_1\\ -q_2& ik& e^{-q_{1}h}{q}_1& -e^{q_{1}h}{q}_1& -e^{-s_{1}h}ik& -i{e}^{s_{1}h}k\\ -k^2{\lambda }_2+{\alpha }_2{q}^2& -2ik{\mu }_2{s}_2& e^{-q_{1}h}(k^2{\lambda }_1-{\alpha }_1{q}_1^2)& -e^{q_{1}h}(-k^2{\lambda }_1+{\alpha }_1{q}_1^2)& 2i{e}^{-s_{1}h}k{\mu }_1{s}_1& -2i{e}^{s_{1}h}k{\mu }_1{s}_1\\ -2ik{\mu }_2{q}_2& -{\mu }_2(k^2+s_2{}^2)& 2i{e}^{-q_{1}h}k{\mu }_1{q}_1& -2i{e}^{q_{1}h}k{\mu }_1{q}_1& e^{-s_{1}h}{\mu }_1(k^2+s_1^2)& e^{s_{1}h}{\mu }_1(k^2+s_1^2)\end{array}\right)$,

 $\det (d)=0$,                                                                                

где ${\alpha }_i={\lambda }_i+2{\mu }_i$, параметры с индексом 1 обозначают величины, относящиеся к слою, а c индексом 2 - к полупространству.

В силу равенства нулю определителя, опустим пятую строку как наиболее громоздкую, а затем перенесем, например, последний столбец в правую часть. Тогда систему уравнений можно представить в следующем виде:

 $\left(\begin{array}{ccccc}0& 0& k^2{\lambda }_1-{\alpha }_1{q}_1^2& -k^2{\lambda }_1+{\alpha }_1{q}_1^2& -2ik{\mu }_1{s}_1\\ 0& 0& -2ik{\mu }_1{q}_1& 2ik{\mu }_1{q}_1& -{\mu }_1(k^2+s_1^2)\\ ik& s_2& -e^{-q_{1}h}ik& -e^{q_{1}h}ik& -e^{-s_{1}h}{s}_1\\ -q{}_2& ik& e^{-q_{1}h}{q}_1& -e^{q_{1}h}{q}_1& -e^{-s_{1}h}ik\\ -2ik{\mu }_2{q}_2& -{\mu }_2(k^2+s_2{}^2)& 2i{e}^{-q_{1}h}k{\mu }_1{q}_1& -2i{e}^{q_{1}h}k{\mu }_1{q}_1& e^{-s_{1}h}{\mu }_1(k^2+s_1^2)\end{array}\right)\left(\begin{array}{c}{A}_2\\ B_2\\ A_1^{(-)}\\ A_1^{(+)}\\ B_1^{(-)}\end{array}\right)=\left(\begin{array}{c}-2ik{\mu }_1{s}_1\\ {\mu }_1(k^2+s_1^2)\\ -e^{s_{1}h}{s}_1\\ i{e}^{s_{1}h}k\\ -e^{s_{1}h}{\mu }_1(k^2+s_1^2)\end{array}\right)B_1^{(+)}$                                                                                                      

Или

$MX=G{B}_1^{(+)}$.

По методу Крамера связь между амплитудами $A_2$ и $B_1^{(+)}$ выражается следующим образом:

$A_2=\frac{\det (M_1)}{\det (M)}{B}_1^{(+)}$,

где

 $M_1=\left(\begin{array}{ccccc}-2ik{\mu }_1{s}_1& 0& k^2{\lambda }_1-{\alpha }_1{q}_1^2& -k^2{\lambda }_1+{\alpha }_1{q}_1^2& -2ik{\mu }_1{s}_1\\ {\mu }_1(k^2+s_1^2)& 0& -2ik{\mu }_1{q}_1& 2ik{\mu }_1{q}_1& -{\mu }_1(k^2+s_1^2)\\ -e^{s_{1}h}{s}_1& s_N& -e^{-q_{1}h}ik& -e^{q_{1}h}ik& -e^{-s_{1}h}{s}_1\\ i{e}^{s_{1}h}k& ik& e^{-q_{1}h}{q}_1& -e^{q_{1}h}{q}_1& -e^{-s_{1}h}ik\\ -e^{s_{1}h}{\mu }_1(k^2+s_1^2)& -{\mu }_N(k^2+s_N{}^2)& 2i{e}^{-q_{1}h}k{\mu }_1{q}_1& -2i{e}^{q_{1}h}k{\mu }_1{q}_1& e^{-s_{1}h}{\mu }_1(k^2+s_1^2)\end{array}\right)$.              

Связь $B_1^{(+)}$ с остальными амплитудами находится аналогичным образом. Итак, амплитуда поверхностной волны оказывается пропорциональна единственной выбранной нормировочной амплитуде (в данном примере $B_1^{(+)}$ ). Поскольку выбор значения нормировочной амплитуды не имеет принципиального значения с точки зрения математических преобразований, то для удобства обозначим ее как F.

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

Приближение плавной горизонтальной неоднородности

Построение поля ПАВ при плавном изменении параметров по горизонтали удобно выполнить в ВКБ-приближении [26]. Представим вектор смещений в следующем виде:

 $u(x,z,t)=V(\epsilon x,z)e^{i\omega (t-kx)},$         (6)

где ε - малый параметр. Для простоты рассмотрим случай однослойного полупространства. Тогда искомое поле удовлетворяет уравнению движения:

 $\begin{array}{l}\rho \frac{{\partial }^{2}u}{\partial t^2}=(\lambda +2\mu )\nabla \text{div}u-\mu \text{rotrot}u+\\ +\nabla \lambda \text{div}u+\nabla \mu \times \text{rot}u+2(\nabla \mu ,\nabla )u\end{array}$    (7)

и граничным условиям:

${\left.T_{zz}\right|}_{z=0}=\mathrm{0,}{\left.T_{xz}\right|}_{z=0}=0$.

Теперь представим амплитудный множитель V в виде ряда по степеням ε:

 $V={\displaystyle \sum _{j=0}^{\infty }{\epsilon }^j{V}^{(j)}}$. (8)

При подставке (6) и (8) в уравнение движения (7) и граничные условия можно убедиться, что в нулевом приближении задача сводится к таковой для горизонтально однородной среды, поэтому необходимо рассмотреть следующий член разложения ε¹. В этом случае для амплитуд компонент смещений может быть получена следующая система уравнений:

 $\begin{array}{c}\frac{\partial }{\partial z}\left((\lambda +2\mu )\frac{\partial V_z^{(1)}}{\partial z}-ik\lambda V_x^{(1)}\right)-ik\mu \frac{\partial V_x^{(1)}}{\partial z}+\\ +\left(\rho {\omega }^2-\mu k^2\right)V_z^{(1)}=L_1\left(V_x^{(0)},V_z^{(0)}\right),\\ \frac{\partial }{\partial z}\left(\mu \frac{\partial V_x^{(1)}}{\partial z}-ik\mu V_z^{(1)}\right)-ik\lambda \frac{\partial V_z^{(1)}}{\partial z}+\\ +\left(\rho {\omega }^2-(\lambda +2\mu )k^2\right)V_x^{(1)}=L_2\left(V_x^{(0)},V_z^{(0)}\right),\end{array}$       (9)

 $\begin{array}{l}\text{(}\lambda +2\mu )\frac{\partial V_z^{(1)}}{\partial z}-ik\lambda V_x^{(1)}=-\frac{\lambda }{c}\frac{\partial V_x^{(0)}}{\partial x},\\ \mu \left(\frac{\partial V_x^{(1)}}{\partial z}-ik{V}_z^{(1)}\right)=-\frac{\mu }{c}\frac{\partial V_z^{(0)}}{\partial x},\end{array}$  (10)

где операторы $L_1\left(V_x^{(0)},V_z^{(0)}\right)$ и $L_2\left(V_x^{(0)},V_z^{(0)}\right)$ имеют вид:

$\begin{array}{l}{L}_1\left(V_x^{(0)},V_z^{(0)}\right)=-\frac{1}{c}\left\{\frac{\partial }{\partial z}\left[\lambda \frac{\partial V_x^{(0)}}{\partial x}\right]+\frac{\partial }{\partial x}\left[\mu \frac{\partial V_x^{(0)}}{\partial z}\right]\right\}+\frac{i}{c}\left\{V_z^{(0)}\frac{\partial }{\partial x}(k\mu )+2\mu k\frac{\partial V_z^{(0)}}{\partial x}\right\},\\ L_2\left(V_x^{(0)},V_z^{(0)}\right)=-\frac{1}{c}\left\{\frac{\partial }{\partial x}\left[\lambda \frac{\partial V_z^{(0)}}{\partial z}\right]+\frac{\partial }{\partial z}\left[\mu \frac{\partial V_z^{(0)}}{\partial x}\right]\right\}+\frac{i}{c}\left\{2(\lambda +2\mu )k\frac{\partial V_x^{(0)}}{\partial x}+V_x^{(0)}\frac{\partial }{\partial x}((\lambda +2\mu )k)\right\}.\end{array}$                                                                                                      

В итоге получена неоднородная краевая задача. Для получения условия ее разрешимости, умножим первое из уравнений (9) на $V_z^{(0)}$, а второе на $V_x^{(0)}{}^{*}$ (звездочка - комплексное сопряжение), вычтем первое из второго и проинтегрируем по z с учетом граничных условий.

Получаем:

 $\frac{\partial }{\partial x}\left[k{I}_1+I_2\right]=\mathrm{0,}$                                                                  (11)

где

$\begin{array}{l}{I}_1={\displaystyle \underset{0}{\overset{\infty }{\int }}\left[\text{(}\lambda +2\mu ){\left(V_x^{(0)}\right)}^2+\mu {\left(V_z^{(0)}\right)}^2\right]}dz,\\ I_2={\displaystyle \underset{0}{\overset{\infty }{\int }}\left[\mu V_z^{(0)}\frac{\partial V_x^{(0)}}{\partial z}-\lambda V_x^{(0)}\frac{\partial V_z^{(0)}}{\partial z}\right]}dz.\end{array}$

Уравнение (11) с точностью до переобозначений можно переписать в следующем виде:

 $\frac{\partial }{\partial x}\left[\omega {ñ}_g{\displaystyle \underset{0}{\overset{\infty }{\int }}\rho \left(V_x^2+V_z^2\right)dz}\right]=\mathrm{0,}$                                            (12)

где $c_g$ - групповая скорость волны.

Величина в квадратных скобках в (12) является потоком энергии через полубесконечное сечение и, как нетрудно показать, эта величина в первом приближении остается постоянной [21]:

 $P={\omega }^2{c}_g(\omega ){\displaystyle \underset{0}{\overset{\infty }{\int }}\rho (z)[V_x^2+V_z^2]}dz=\text{const}$.           (13)

Отметим, что соотношение (13) не является впервые полученным в настоящее работе. Данные выкладки приведены для более полного описания используемой математической модели.

Обобщение соотношения (13) для многослойной системы можно представить следующим образом:

$P={\omega }^2{c}_g(\omega ){\displaystyle \sum _{i=1}^N{\displaystyle \underset{0}{\overset{h_i}{\int }}{\rho }_i[V_{xi}^2+V_{zi}^2]}dz}$.

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

Так как величина потока энергии пропорциональна сумме квадратов компонент смещений, а они в свою очередь пропорциональны единственной выбранной нормировочной амплитуде, то поток энергии пропорционален квадрату этой нормировочной амплитуды: P ~ F².

Для случая однородного полупространства (однослойной среды) выражения для компонент смещений можно выразить в явном виде [22]:

 $\begin{array}{l}{u}_x=Fk(e^{-qz}-\frac{2qs}{k^2+s^2}{e}^{-sz})e^{ikx},\\ u_z=iFq(e^{-qz}-\frac{2k^2}{k^2+s^2}{e}^{-sz})e^{ikx}.\end{array}$             (14)

Тогда для потока энергии (13) имеем:

$\begin{array}{c}P(c_l,c_t,\rho )={\omega }^2{c}_g\rho \left(\frac{k^2-q^2}{2q}-\frac{4k^2{q}^2(k^2-s^2)}{2s{\left(k^2+s^2\right)}^2}\right)F^2=\\ =C(\omega ,c_l,c_t,\rho )F^2.\end{array}$

Рассмотрим теперь две среды, представляющие собой полупространства с постоянными значениями упругих параметров. Тогда для первой среды $P_1(\omega ,c_l,c_t,\rho )=C_1(\omega ,c_l,c_t,\rho )F_1{}^2$, а для второй $P_2(\omega ,c_l,c_t,\rho )=C_2(\omega ,c_l,c_t,\rho )F_2{}^2$. Будем считать, что эти полупространства связаны областью с малыми градиентами упругих параметров, обеспечивающих плавный переход ПАВ из одной среды в другую без отражений, если $P_1(\omega ,c_l,c_t,\rho )=P_2(\omega ,c_l,c_t,\rho )$. Отсюда следует, что отношение нормировочной амплитуды ПАВ во второй среде к амплитуде в первой среде выражается следующим образом

$\frac{F_2}{F_1}=\sqrt{\frac{C_1(\omega ,c_l,c_t,\rho )}{C_2(\omega ,c_l,c_t,\rho )}}$.    (15)

Из этого выражения следует, что, зная частоту и параметры двух сред, можно определить отношение нормировочных амплитуд $F_2/F_1$. Поскольку амплитуда компонент смещений в конкретной среде пропорциональна соответствующей нормировочной амплитуде (14), то относительное изменение будет пропорционально $F_2/F_1$. Подставляя (15) и параметры сред в (14), можно получить искомые изменения амплитуды компонент смещений.

При рассмотрении плоской по горизонтали поверхностной волны можно получить только соотношения между амплитудами, но не абсолютные значения, которые возможно отыскать только при добавлении источника возмущения, например, точечного [27]. Однако в этом случае необходимо учитывать, во-первых, влияние объемных волн и, во-вторых, тот факт, что на поверхностную волну идет не вся энергия источника, а только ее часть, причем в зависимости от параметров системы, эта доля может меняться [20]. Это усложняет анализ амплитудных закономерностей поверхностной волны. Поэтому, на наш взгляд, представляется более целесообразным развивать описание в отсутствие источника, поскольку во многих практических случаях интерес представляет именно относительное изменение, а не абсолютное значение.

Численная модель

Для проверки возможности подобной аналитической оценки амплитуды ПАВ использовалось численное моделирование методом конечных элементов в программном комплексе COMSOL Multiphysics 5.3 с подключенным модулем Solid Mechanics (лицензия № 9600341). Главные отличия численного моделирования от аналитического рассмотрения заключаются в необходимости использовать источник сигнала и рассматривать модель ограниченного размера, который, однако, может быть довольно велик по сравнению с длиной используемой волны, поскольку расчет в двумерной постановке в пространстве частот не слишком требователен к вычислительным ресурсам и может быть выполнен на современных персональных компьютерах за минуты или часы.

Основные использовавшиеся принципы построения численной модели представлены в [28, 29]. Геометрия модели рассматриваемой среды представляет собой вытянутый вдоль оси x прямоугольник, разбитый на горизонтальные слои, которые могут отличаться между собой упругими параметрами, а также состоящий из нескольких вертикальных зон, часть из которых представляют собой плоскослоистые участки среды, рассмотренные в предыдущем разделе, а другие - технические зоны, обеспечивающие плавный переход упругих параметров от различных плоскослоистых участков, а также стартовую и финишную зоны, необходимые для обеспечения генерации и поглощения ПАВ (рис. 2).

Рис. 2. Пример модели среды, сетки конечных элементов, использовавшихся при численном моделировании, и результатов расчета распределения амплитуды ПАВ вдоль свободной границы среды.

Упругие параметры (плотность и скорости продольной и поперечной волн) в стартовой и финишной зонах не меняются с глубиной, несмотря на то, что присутствует формальное разделение на слои (подобное деление удобно для автоматического построения сетки конечных элементов). В этом случае в качестве источника ПАВ мы можем использовать существующее [22] аналитическое решение для волны Рэлея с известной длиной, которая используется в качестве характерного линейного размера модели (рис. 2).

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

Для трансформации ПАВ без отражений между зонами (сначала переход волны Рэлея из однородного полупространства в зоне 0 в слоистую среду в зоне 1, затем - в другую слоистую среду в зоне 2, а после - в финишную однородную зону) используются протяженные зоны с плавным изменением упругих параметров от значений в левой зоне к значению в правой (для каждого слоя отдельно).

Протяженность каждой из зон 0-3 составляет 30 длин ПАВ для стартовой зоны и достаточна для имитации бесконечной по x среды. Зоны плавного перехода упругих параметров шире (50 длин ПАВ для стартовой зоны), что обеспечивает трансформацию ПАВ без видимых отражений. Глубина моделируемой среды ограничена семью длинами ПАВ для стартовой зоны. Отметим, что эти параметры (как и сетка конечных элементов) избыточны и позволяют получить решения с ошибкой не более 0.1%. Эта оценка была получена стандартной проверкой сходимости численного решения путем серии последовательных расчетов при сгущении сетки конечных элементов и увеличении линейных размеров зон моделируемой среды. Валидация путем сравнения результата расчета характеристик волны Рэлея для случая однородной среды с известным аналитическим решением доказала корректность составленной численной модели, поскольку полученные решения совпали с точностью до четвертого знака после запятой.

Физические процессы в численном моделировании также как и в теоретическом рассмотрении описываются уравнениями (1) механики идеальной линейной упругой среды. При этом граничное условие при z = 0 соответствует свободной поверхности, на нижнюю и правую границу наложены слабоотражающие условия, а на левой границе при x = 0 расположен источник, задаваемый смещениями (14).

Для большего ослабления эффектов от ограниченного размера численной модели вдоль правой и нижней границ расположены идеально согласованные слои (Perfectly Matched Layers - PMLs), поглощающие большую часть падающих на них упругих волн.

Разбиение на конечные элементы (рис. 2) выполнено с помощью структурированной сетки (все элементы представляют собой прямоугольники) с линейным размером элементов, не превышающим 0.1 длины ПАВ для стартовой зоны. Отметим, что упругие параметры стартовой зоны выбираются таким образом, чтобы скорость ПАВ, а соответственно, и ее длина были наименьшими из возможных в рассматриваемой системе, благодаря этому длины всех распространяющихся в данной системе упругих волн превосходят размер конечных элементов не менее, чем в 10 раз. Вблизи свободной границы модели для еще большего повышения точности решения происходит сгущение сетки конечных элементов с уменьшением вертикального размера конечного элемента до 5% от первоначального размера в геометрической прогрессии с коэффициентом 1.2.

Решение численной модели происходит в режиме варьирования параметров (автоматический последовательный расчет нескольких численных реализаций, отличающихся значениями одного или нескольких входных параметров модели) с помощью прямого решателя MUMPS (MUltifrontal Massively Parallel Solver) [30], реализующего LU-факторизацию.

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

Результаты

Перейдем к рассмотрению полученных результатов. Начнем с наиболее простой модели - перехода из однослойной системы в однослойную, т.е. трансформации волны Рэлея. В дальнейшем первою среду (из которой происходит переход) будем называть опорной, а вторую среду (в которую происходит переход) будем называть исследуемой. При этом будем рассматривать зависимости изменения амплитуды ПАВ на поверхности (поскольку именно эта величина как правило измеряется в эксперименте) от коэффициента изменения K, т.е. величины, отражающей во сколько раз значение того или иного параметра исследуемой среды отличается от соответствующего значения опорной среды: $K_{\rho }$ - коэффициент изменения плотности, $K_{c_l}$ - коэффициент изменения скорости продольных волн, $K_{c_t}$ - коэффициент изменения скорости поперечных волн. При варьировании какого-либо коэффициента K остальные коэффициенты принимаются равными 1; например, при изменении $K_{\rho }$ считаем, что $K_{c_t}=K_{c_l}=1$. Параметры рассматриваемых сред представлены в табл. 1.

Таблица 1. Параметры опорной и исследуемой среды.

На рис. 3 представлены зависимости изменения амплитуды горизонтальной и вертикальной компоненты смещений на поверхности от коэффициента изменения различных упругих параметров, полученные аналитически и с помощью численного моделирования. При увеличении коэффициента изменения плотности и скорости продольных волн амплитуда ПАВ уменьшается. В случае изменения скорости поперечных волн наблюдается локальный минимум. Это связано с тем, что с увеличением скорости поперечных волн увеличивается роль векторного потенциала [31] и, как следствие, увеличение амплитуды на поверхности. Кроме того, отметим, что для однослойной модели наиболее сильное изменение амплитуды связано с изменением скорости продольных волн. В целом разница теории и численного моделирования составляет менее 1.5%, однако в области низких с_l и высоких с_t начинается расхождение, связанное с уменьшением коэффициента Пуассона. Однако поскольку для большинства сплавов, металлов, горных пород значение коэффициента Пуассона лежит в пределах 0.25-0.35 [32], то области вне этих пределов выходят за рамки данного исследования.

Рис. 3. Зависимость относительного изменения амплитуды горизонтальной и вертикальной компонент смещений на поверхности от (а) – коэффициента изменения плотности, (б) – скорости продольных волн, (в) – скорости поперечных волн, (г) – совместного изменения скорости продольных и поперечных волн, (д) – совместного изменения всех трех упругих параметров; сплошная линия – теория, черные точки – численное моделирование (зона перехода – 50λ), серые точки – численное моделирование (зона перехода – λ).

При совместном изменении скоростей или всех трех упругих параметров, разница теоретических и численных расчетов становится менее 0.5%. Кроме того, отметим, что при переходе в более “мягкую” среду наблюдается увеличение амплитуды, а при переходе в более “твердую” - уменьшение, что согласуется с основами ММЗ [3].

Перейдем к более сложной модели - переходу из однослойной среды в двухслойную. Выберем в качестве исходной двухслойной модели скоростную модель PREM [33] с усредненными по 24.5 км значениями, а для однослойной модели выберем параметры верхнего слоя, приведенные в табл. 2.

Таблица 2. Параметры опорной и исследуемой среды.

Введем следующие обозначения: $K_{c_l{}_1}=K_{c_t{}_1}\equiv K_{c_1}$ - коэффициент совместного изменения скорости продольных и поперечных волн в верхнем слое относительно исходной модели, $K_{c_l{}_2}=K_{c_t{}_2}\equiv K_{c_2}$ - коэффициент совместного изменения скорости продольных и поперечных волн в нижнем слое относительно исходной модели, $K_H$ - коэффициент изменения толщины слоя относительно исходной модели, $H=24.5\text{}\text{êì}$.

На рис. 4 представлены зависимости изменения амплитуды при совместном изменении скорости продольных и поперечных волн в первом и втором слое, а также при изменении толщины слоя. При увеличении скоростей в верхнем слое амплитуда уменьшается, а при изменении скорости в нижнем слое увеличивается. В зависимости от мощности слоя наблюдается максимум. Отметим, что в предельных случаях большой и малой мощности двухслойная система эквивалентна однослойной, поэтому зависимости асимптотически стремятся к постоянной величине. В данном случае при увеличении K_H эта величина равна 1, поскольку свойства опорной среды такие же, как и у верхнего слоя исследуемой среды. Кроме того, в области малых толщин наблюдается более заметное расхождение теории с численным моделированием. Предположительно, это связано с тем, что в этой области существует дополнительная неоднородная волна, распространяющаяся вдоль границы слоя и полупространства [22], которая учитывается в численном моделировании, но игнорируется в теоретических расчетах.

Рис. 4. Зависимость относительного изменения амплитуды горизонтальной и вертикальной компонент смещений на поверхности от коэффициента совместного изменения скорости продольных и поперечных волн (а) – в верхнем и (б) – нижнем слое, (в) - коэффициента изменения мощности слоя; λ/H = 2.5; сплошная линия – теория, черные точки – численное моделирование (зона перехода – 50λ), серые точки – численное моделирование (зона перехода – λ).

Отметим, что из рис. 3 и 4 следует, что при замене плавного перехода размером 50 длин волн на резкий переход (что более характерно для реальных ситуаций) размером 1 длина волны наблюдается большее расхождение с теорией. Это связано с тем, что в данном случае часть энергии отражается от границы [34, 35] и условие сохранения потока энергии не выполняется. Однако оценки показывают, что для рассматриваемых случаев это расхождение оказывается около 5%, поэтому можно считать, что при небольших контрастах результаты расчетов для плавного перехода можно применять и для резкого.

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

Рис. 5 демонстрирует, что при определенных длинах волн наблюдается максимум изменения амплитуды, минимумы групповой скорости и эллиптичности. При этом их положения зависят как от толщины слоя, так и от упругих свойств среды. Для рассматриваемых случаев они расположены в пределах от 1.5Н до 2.5Н. Это согласуется с ММЗ, согласно которому наибольшая чувствительность соответствует глубинам неоднородности, равным примерно 0.3-0.5 длины волны [8].

Рис. 5. Зависимость (а) – относительного изменения амплитуды горизонтальной и вертикальной компонент смещений на поверхности, (б) – фазовой и групповой скорости и (в) – эллиптичности на поверхности от длины волны; сплошная черная – K_c1 = 0.8, K_c2 = 1, K_H = 1, сплошная серая – K_c1 = 1, K_c2 = 0.8, K_H = 1, штриховая черная – K_c1 = 1.2, K_c2 = 1, K_H = 1, штриховая серая – K_c1 = 1, K_c2 = 1.2, K_H = 1, штрихпунктирная черная – K_c1 = 1, K_c2 = 1, K_H = 1.2, штрихпунктирная серая – K_c1 = 1, K_c2 = 1, K_H = 0.8.

На рис. 6 представлены зависимости доминантной длины волны ${\lambda }_0$, т.е. длины волны, при которой наблюдается экстремум (см. рис. 5) от различных коэффициентов изменения параметров среды. При увеличении скорости продольных и поперечных волн в нижнем слое, а также мощности, наблюдается увеличение доминантной длины волны. В случае же увеличения скорости продольных и поперечных волн в верхнем слое наблюдается уменьшение. Отсюда следует, что доминантная длина волны тем больше, чем выше контраст между слоями и чем больше толщина слоя. В работе [8] были произведены расчеты коэффициента глубинной привязки ММЗ, т.е. величины, аналогичной доминантной длине волны, от скоростного контраста неоднородности и было получено, что эта величина возрастает с увеличением контраста, что согласуется с рис. 6.

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

Перейдем к рассмотрению перехода ПАВ из двухслойной среды в двухслойную. Зафиксируем параметры опорной среды (табл. 3) и, как и в предыдущем случае, будем варьировать параметры исследуемой среды: верхнего и нижнего слоя, а также мощность верхнего слоя. Очевидно, что зависимости от изменения упругих параметров будут аналогичны рис. 4, поскольку исследуемая среда не изменилась. Поэтому ограничимся только зависимостями от длины волны.

Таблица 3. Параметры опорной и исследуемой среды.

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

Рис. 7. Зависимость относительного изменения амплитуды горизонтальной и вертикальной компонент смещений на поверхности от длины волны; сплошная черная – K_c1 = 0.8, K_c2 = 1, K_H = 1, сплошная серая – K_c1 = 1, K_c2 = 0.8, K_H = 1, штриховая черная – K_c1 = 1.2, K_c2 = 1, K_H = 1, штриховая серая – K_c1 = 1, K_c2 = 1.2, K_H = 1, штрихпунктирная черная – K_c1 = 1, K_c2 = 1, K_H = 1.2, штрихпунктирная серая – K_c1 = 1, K_c2 = 1, K_H = 0.8.

Заключение

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

Применительно к задачам пассивной поверхностно-волновой сейсморазведки результаты данной работы позволяют сделать следующие выводы. При переходе из однослойной системы в однослойную зависимости амплитуды от изменения плотности и скорости продольных волн соответствуют ожиданиям, в то время как зависимость от скорости поперечных волн имеет иной вид, что говорит о более сложной зависимости свойств поверхностной волны от этого параметра. Изменение амплитуды в этом случае также меньше, чем в зависимости от скорости продольных волн. Это важный результат, поскольку на данный момент в ММЗ полагается, что коэффициент Пуассона постоянен, т.е. не выделяется, какая из скоростей упругих волн вносит определяющий вклад, но основное внимание уделяется скорости поперечных волн [36]. Однако было показано, что, в частности, для однородного полупространства вклад скорости продольных волн существенно больше, поэтому необходимо дальнейшее изучение этого вопроса. При переходе в двухслойную среду зависимости от упругих параметров, а также частотные зависимости демонстрируют соответствие ММЗ. Выбор опорной среды влияет на значение доминантной длины волны, которая может меняться в пределах от 0.3 до 0.7 мощностей слоя, что связано с изменением профиля ПАВ.

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

Работа поддержана грантом Российского научного фонда № 22-77-00067.