Tidal generated electric field in the multi-layer structure and the possibilities of its employment for deriving the elastic properties and permeability of the subsurface formations

Full Text

1. ВВЕДЕНИЕ

Электрокинетический механизм возникновения электрического поля в горных породах, испытывающих упругие деформации, описывается в рамках уравнения Гельмгольца$–$Смолуховского в сочетании с теорией Био пороупругой флюидонасыщенной среды [Френкель, 1944; Biot, 1956; Светов, 2008]. В качестве источника таких деформаций могут выступать волны лунно-солнечных приливов, вызывающих периодические колебания порового давления и, как следствие, электрокинетических потенциалов. Соответствующая теория для этого случая приводится в работах [Garagash et al., 2005; Гохберг и др., 2007; Гохберг и др., 2009], где получено аналитическое решение задачи для распределения порового давления в однородной и двухслойной сферически-симметричных глобальных моделях, а также плоской однородной модели среды. В качестве важной особенности этого распределения отмечается наличие скин-слоя вблизи границ сред с различными упругими и петрофизическими параметрами [Гохберг и др., 2009]. Установлено, что в латерально-однородных средах информация о проницаемости и упругих параметрах содержится только в вертикальной компоненте электрического поля [Гохберг и др., 2007].

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

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

1.1. Обобщение решения в сферической модели на случай n слоев

Рассмотрим глобальную сферически-симметричную модель Земли, состоящую из n сферических слоев, границы между которыми определяются набором значений r₀, r₁, r₂,..., r_n$–$1, r_n, n > 1 где r₀ = a, r_n = 0. Воспользуемся полным представлением для порового давления в слое n сферической модели и ограничимся рассмотрением полусуточной приливной гармоники M2 [Гохберг и др., 2009]:

 

$$\begin{array}{c}{p}_n(r,\theta ,\lambda )={\displaystyle \sum _{m=1}^{\infty }{D}_m^{(n)}{j}_m(z){\displaystyle \sum _{k=-m}^m{C}_{mk}^{(n)}}{Y}_m^{(k)}}-\\ -\frac{q_n}{d_n}G{\left(\frac{r}{a}\right)}^2{\sin }^2\theta e^{i2\lambda }.\end{array}$$

(1)

Здесь: p_n $–$ поровое давление в слое n, представленное в сферической системе координат $r,\theta ,\lambda ;$ $j_m(z)$ $–$ функции Бесселя; $Y_m^{(k)}$ $–$ сферические функции.

Величины $\frac{q_n}{d_n}G$ в выражении (1) зависят от упругих свойств слоев и объемной деформации слоев (константы $\frac{q_n}{d_n}$ определяются в работе [Гохберг и др., 2009]; G $–$ постоянная Дудсона из теории приливов [Сидоренков, 2002]). В случае полусуточной приливной волны M2 с частотой $\omega \approx \frac{2\pi }{12}$ рад/ч, все постоянные коэффициенты $C_{mk}^{\left(n\right)}$ с m ≠ 2 и k ≠ ±2, а также все коэффициенты $D_m^{\left(n\right)}$ с m ≠ 2 должны обращаться в нуль, чтобы выполнялось условие равенства давления нулю на земной поверхности. В результате, выражение (1) приводится к виду:

$$\begin{array}{c}{p}_n(r,\theta ,\lambda )=\\ =\left\{D_2^{(n)}{C}_{2\pm 2}^{(n)}{j}_2(r/{\delta }_n\text{})-\frac{q_n}{d_n}G{\left(\frac{r}{a}\right)}^2\right\}{\sin }^2\theta e^{i2\lambda }.\end{array}$$

(2)

Здесь величина ${\delta }_n=\sqrt{\frac{2{\kappa }_n{M}_n}{{\eta }_n\omega }},$ зависящая от вязкоcти поровой жидкости η, коэффициента проницаемости k и упругих параметров среды α, M (коэффициента Био и модуля Био), а также частоты приливной волны ω, и имеющая размерность расстояния, определяет обратный скин-эффект. Аргументом функций Бесселя является величина (1$–$i)r/δ, однако для краткости в последующем изложении мы опускаем множитель (1$–$i).

На границах слоев выполняется условие непрерывности порового давления:

pn = pn$–$1 при r = rn$–$1, при n > 1; p0 = 0 при r = r0. (3)

Произведение коэффициентов $D_2^{\left(n\right)}{C}_{2\pm 2}^{\left(n\right)},$ определяется из условия непрерывности давления на границах слоев путем прогонки, начиная с верхней границы r0:

Cлой 1: $D_2^{(1)}\text{}\cdot C_{2\pm 2}^{(1)}=\frac{q_1}{d_1}G\frac{1}{j_2(a/{\delta }_1)}.$

Cлой 2: $D_2^{(2)}\cdot C_{2\pm 2}^{(2)}=\frac{q_1}{d_1}G\frac{j_2(r_1\text{}/\text{}{\delta }_1)}{j_2(r_0/{\delta }_1\text{})j_2\text{}(r_1/{\delta }_2\text{})}+$

$+\left(\frac{q_2}{d_2}-\frac{q_1}{d_1}\right)G\frac{1}{j_2(r_1/{\delta }_2)}{\left(\frac{r_1}{a}\right)}^2.$

Cлой 3:

$D_2^{(3)}\cdot C_{2\pm 2}^{(3)}=G\left\{\frac{q_1}{d_1\text{}}\frac{j_2(r_1/{\delta }_1)j_2(r_2/{\delta }_2\text{})}{j_2(r_0/{\delta }_1\text{})j_2(r_1/{\delta }_2\text{})j_2(r_2/{\delta }_3\text{})}+\right.$

$\begin{array}{l}+\left(\frac{q_2}{d_2}-\frac{q_1}{d_1}\right)\frac{j_2(r_2/{\delta }_2)}{j_2(r_1/{\delta }_2)j_2(r_2/{\delta }_3)}{\left(\frac{r_1}{a}\right)}^2+\\ \left.+\left(\frac{q_3}{d_3}-\frac{q_2}{d_2}\right)\frac{1}{j_2(r_2/{\delta }_3)}{\left(\frac{r_2}{a}\right)}^2\right\}.\end{array}$

Слой n:

$D_2^{(n)}\cdot C_{2\pm 2}^{(n)}=G\left(\frac{q_1}{d_1}\frac{{\displaystyle \prod _{k=1}^{n-1}{j}_2(r_k/{\delta }_k)}}{{\displaystyle \prod _{k=1}^n{j}_2(r_{k-1}/{\delta }_k)}}+\right.$

 

$+\left.{\displaystyle \sum _{l=1}^{n-1}\left(\frac{q_{l+1}}{d_{l+1}}-\frac{q_l}{d_l}\right)\frac{{\displaystyle \prod _{k=l+1}^{n-1}{j}_2(r_k/{\delta }_k)}}{{\displaystyle \prod _{k=l+1}^n{j}_2(r_{k-1}/{\delta }_k)}}{\left(\frac{r_l}{a}\right)}^2}\right).$ (4)

При этом, мы определяем операцию произведения таким образом, что:

${\displaystyle \prod _{i=a}^b{x}_i}=\mathrm{1,}$ если a > b.

Таким образом, решение в слое n имеет вид:

$p_n(r,\theta ,\lambda )=G\left\{\left(\frac{q_1}{d_1}\frac{{\displaystyle \prod _{k=1}^{n-1}{j}_2(r_k/{\delta }_k)}}{{\displaystyle \prod _{k=1}^n{j}_2(r_{k-1}/{\delta }_k)}}\right.\right.\left.+{\displaystyle \sum _{l=1}^{n-1}\left(\frac{q_{l+1}}{d_{l+1}}-\frac{q_l}{d_l}\right)\frac{{\displaystyle \prod _{k=l+1}^{n-1}{j}_2(r_k/{\delta }_k)}}{{\displaystyle \prod _{k=l+1}^n{j}_2(r_{k-1}\text{}/\text{}{\delta }_k)}}{\left(\frac{r_l}{a}\right)}^2}\right)\text{}{j}_2(r/{\delta }_n\text{})-\frac{q_n}{d_n}{\left(\frac{r}{a}\right)}^2\left.\underset{}{\overset{}{\underset{}{\overset{}{\underset{}{\overset{}{{}_{}^{}}}}}}}\right\}{\sin }^2\text{}\theta e^{i2\lambda }.$ (5)

Проведем некоторые преобразования (5), чтобы представить его в более простом виде. Воспользуемся асимптотикой функций Бесселя, справедливой при больших значениях r/δ:

$j_2\left((1-i)\raisebox{1ex}{$r$}\!\left/ \!\raisebox{-1ex}{$\delta $}\right.\right)\cong \frac{i}{2(1-i)\raisebox{1ex}{$r$}\!\left/ \!\raisebox{-1ex}{$\delta $}\right.}\exp \left[(1+i)\raisebox{1ex}{$r$}\!\left/ \!\raisebox{-1ex}{$\delta $}\right.\right].$

Тогда можно видеть, что

$$p_n(r,\theta ,\lambda )=G{\sin }^2\theta e^{i2\lambda }\left\{\left.{\displaystyle \sum _{l=0}^{n-1}\left(\left(\frac{q_{l+1}}{d_{l+1}}-\frac{q_l}{d_l}\right)\exp \left\{\left(1+i\right)\left({\displaystyle \sum _{k=l+1}^n\frac{r_k-r_{k-1}}{{\delta }_k}}+\frac{r}{{\delta }_n}\right)\right\}{\left(\frac{r_l}{a}\right)}^2\right)}\right)-\frac{q_n}{d_n}{\left(\frac{r}{a}\right)}^2\right\}.$$

(6)

При записи выражения (6) мы воспользовались обозначением $\frac{q_0}{d_0}=0.$

Отметим, что для общности записи мы использовали символы суммирования, которые имеют смысл лишь при l ≤ n$–$1, k ≤ n, а в противном случае соответствующие слагаемые в выражениях отсутствуют.

Переобозначив константы, характеризующие свойства слоев и деформацию:

$$\frac{q_n}{d_n}G=\frac{1}{2}{\alpha }_n{M}_n{\epsilon }_0,$$

(7)

и, считая, что α0M0 = 0, запишем решение:

$p_n(r,\theta ,\lambda )=\frac{1}{2}{\epsilon }_0{\sin }^2\theta e^{i2\lambda }\text{}\left\{{\displaystyle \sum _{l=0}^{n-1}\left(\left({\alpha }_{l+1}{M}_{l+1}\text{}-{\alpha }_l{M}_l\right)\exp \left\{\left(1+i\right)\left({\displaystyle \sum _{k=l+1}^n\frac{r_k-r_{k-1}}{{\delta }_k}}+\frac{r}{{\delta }_n}\right)\right\}{\left(\frac{r_l}{a}\right)}^2\right)-}{\alpha }_n{M}_n{\left(\frac{r}{a}\right)}^2\right\}$

(8)

Дальнейший анализ будем проводить в терминах модуля Био M, коэффициента Био α и коэффициента проницаемости k.

1.2. Переход к горизонтально-слоистой модели

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

Зависимость приливной деформации от широты очень проста, и в случае гармоники M₂ сводится к наличию в выражении для порового давления множителя sin2θ, где θ $–$ коширота (полярный угол, отсчитываемый в меридиональной плоскости) (см. формулу (2)). Для определенности всюду ниже мы ограничиваемся случаем θ = π/2, то есть, рассматривается поле на экваторе. Очевидно, что при работе с реальными данными учет широты не будет представлять сложности. Приливный потенциал для волны O1 также характеризуется несложной зависимостью от кошироты, что потребует лишь незначительной модификации соответствующих выражений.

Таким образом, переходя от сферической модели к локально-плоской модели с учетом соотношений $\frac{r_l}{a}\approx \mathrm{1,}\frac{r}{a}\approx \mathrm{1,}$ $e^{i2\lambda }=\exp [i(\omega t+\phi )],$ исключая из рассмотрения широту и вводя глубину z, отсчитываемую от поверхности земли: $z=a-r,$ а также, имея ввиду, что

 

$$-{\displaystyle \sum _{k=l+1}^n\frac{h_k}{{\delta }_k}}+\frac{a-z}{{\delta }_n}=-{\displaystyle \sum _{k=l+1}^{n-1}\frac{h_k}{{\delta }_k}}-\frac{1}{{\delta }_n}\left(z-{\displaystyle \sum _{k=1}^{n-1}{h}_k}\right),$$

(9)

получаем:

 

$$p_n(z)=\frac{1}{2}{\epsilon }_0{e}^{i\omega t}\left\{{\displaystyle \sum _{l=0}^{n-1}\left(\left({\alpha }_{l+1}\text{}{M}_{l+1}\text{}-{\alpha }_l{M}_l\right)\exp \left\{\left(1+i\right)\left(-{\displaystyle \sum _{k=l+1}^{n-1}\frac{h_k}{{\delta }_k}}-\frac{1}{{\delta }_n}\left(z\text{}-{\displaystyle \sum _{k=1}^{n-1}{h}_k}\right)\right)\right\}\right)}-{\alpha }_n{M}_n\right\}.$$

(10)

Переходя к действительным величинам (во временную область), будем иметь:

 

$p_n(z,t)=\frac{1}{2}{\epsilon }_0\left\{{\displaystyle \sum _{l=0}^{n-1}\left({\alpha }_{l+1}{M}_{l+1}\text{}-{\alpha }_l{M}_l\right)\exp \left(-{\displaystyle \sum _{k=l+1}^n\frac{h_k}{{\delta }_k}+\frac{a-z}{{\delta }_n}}\right)\cos \left(\omega t-\text{}{\displaystyle \sum _{k=l+1}^n\frac{h_k}{{\delta }_k}+\frac{a-z}{{\delta }_n}}\right)}-\text{}{\alpha }_n{M}_n\text{}\cos \left(\omega t\right)\right\},$ (11)

где $h_k=r_{k-1}-r_k$ $–$ мощность слоя k.

Например, для расчета порового давления в первом и втором (сверху) слоях можно использовать следующие выражения:

 

$p_1(z,t)=\frac{1}{2}{\alpha }_1{M}_1{\epsilon }_0\left[\exp \left(-\frac{z}{{\delta }_1}\right)\cos \left(\frac{z}{{\delta }_1}-\omega t\right)-\cos \left(\omega t\right)\right],$  (12)

 

$$\begin{array}{c}{p}_2(z,t)=\frac{1}{2}{\alpha }_2{M}_1{\epsilon }_0\left\{e^{-(z/{\delta }_2-h/{\delta }_2+h/{\delta }_2)}cos(\omega t-z/{\delta }_2+h/{\delta }_2-h/{\delta }_1)\right.+\\ +\left(\frac{{\alpha }_2{M}_2}{{\alpha }_1{M}_1}-1\right)e^{-(z/{\delta }_2-h/{\delta }_2)}\cos (\omega t-z\text{}/\text{}{\delta }_2+h/{\delta }_2)-\frac{{\alpha }_2{M}_2}{{\alpha }_1{M}_1}\cos (\omega t)\left.{}_{{}_{}^{}}^{}\right\}.\end{array}$$

(13)

В целях более компактной записи, дальнейший анализ мы будем вести в терминах комплексных амплитуд (в частотной области).

Дифференцируя (10) по z, получаем выражение для вертикальной составляющей градиента порового давления:

 

$$\frac{\partial p_n(z)}{\partial z}=-\frac{\left(1+i\right)}{2{\delta }_n}{\epsilon }_0{e}^{i\omega t}{\displaystyle \sum _{l=0}^{n-1}\left(\left({\alpha }_{l+1}{M}_{l+1}-{\alpha }_l{M}_l\text{}\right)\exp \left\{\left(1+i\right)\left(-{\displaystyle \sum _{k=l+1}^{n-1}\frac{h_k}{{\delta }_k}}-\frac{1}{{\delta }_n}\left(z-{\displaystyle \sum _{k=1}^{n-1}{h}_k}\right)\right)\right\}\right)}$$

(14)

Из (10) и (14) видно, что решение в каждом последующем слое (от поверхности вниз) есть сумма решений для полупространств, параметры которых определяются характеристиками всех вышележащих слоев с соответствующими сдвигами по глубине. Например, решение в первом слое совпадает с решением для полупространства с параметрами ${\alpha }_1,M_1,{\delta }_1,$ а при переходе во второй слой к нему добавляется дополнительное слагаемое. Аналогичным образом, при переходе в нижележащие слои происходит суммирование с соответствующими дополнительными слагаемыми. При этом, амплитуды последних зависят от контраста ${\alpha }_{l+1}{M}_{l+1}-{\alpha }_l{M}_l,$ а также отношений мощностей и скин-толщин вышележащих слоев. Зависимость от глубины в любом слое является экспоненциальной, что делает целесообразным представление амплитуд градиента порового давления в логарифмическом масштабе. В таком представлении $\frac{\partial p(z)}{\partial z}$ имеет вид кусочно-линейной функции глубины z. В качестве примера приведем графики зависимости модуля комплексной амплитуды давления и его вертикального градиента от глубины для ряда трехслойных моделей c различными параметрами слоев (рис. 1).

2. Чувствительность к параметрам горизонтально-слоистой среды

Особенностью полученных выше выражений (10) и (14) является то, что решение в каждом слое не зависит от свойств нижележащей среды. Это обстоятельство не позволяет проводить изучение глубоких слоев по измерениям на поверхности, однако, проводя измерения внутри среды (в скважинах), достаточно легко получить оценки ее геомеханических и петрофизических характеристик, входящих в выражения (10) и (14). В этой связи, представляется целесообразным проанализировать чувствительность градиента порового давления к упругим параметрам и проницаемости слоев.

2.1. Чувствительность к изменению значений коэффициента Био α

Наиболее простой вид имеет выражение для производной вертикального градиента порового давления в слое n по параметру α слое m (m и n могут быть равны):

$\begin{array}{c}\frac{\partial }{\partial {\alpha }_m}\left(\frac{\partial p_n(z)}{\partial z}\right)=-\frac{\left(1+i\right)}{2{\delta }_n}{\epsilon }_0{e}^{i\omega t}\times \\ \times M_m\text{}\exp \left(-\frac{h_m}{{\delta }_m}\text{}-1\right)\exp \left\{\left(1\text{}+i\right)\left(-\text{}\text{}\text{}{\displaystyle \sum _{k=l+1}^{n-1}\frac{h_k}{{\delta }_k}}\text{}-\frac{1}{{\delta }_n}\left(z\text{}\text{}\text{}-\text{}{\displaystyle \sum _{k=1}^{n-1}{h}_k}\right)\right)\right\}.\end{array}$

(15)

Как мы видим, в выражение для чувствительности к коэффициенту Био αm в виде сомножителя входит модуль Био Mm. С учетом широкого диапазона возможных изменений последнего, можно сделать вывод о том, что амплитуда чувствительности к αm может меняться на несколько порядков.

2.2. Чувствительность к изменению значений модуля Био M и коэффициента проницаемости k

Учитывая некоторые различия в функциональной зависимости $\frac{\partial p_n(z)}{\partial z}$ от параметров текущего и вышележащих слоев, при расчете чувствительности целесообразно отдельно рассмотреть два случая: m ≠ n и m = n.

 

Рис. 1. Зависимость амплитуд порового давления (а) и его вертикального градиента (б) от глубины для набора трехслойных моделей. Параметры моделей приводятся в табл. 1. Величина приливной деформации была принята равной ε0 = 2 ∙ 10–8, частота ω = 2π/12 рад/час.

 

Случай 1 (m ≠ n):

$\begin{array}{c}\frac{\partial }{\partial M_m}\left(\frac{\partial p_n(z)}{\partial z}\right)=\\ =-\frac{\left(1+i\right)}{2{\delta }_n}{\epsilon }_0{e}^{i\omega t}{\displaystyle \sum _{l=0}^{n-1}\exp \left\{\left(1+i\right)\text{}\left(-{\displaystyle \sum _{k=l+1}^{n-1}\frac{h_k}{{\delta }_k}\text{}\text{}}-\text{}\frac{1}{{\delta }_n}\left(z-\text{}{\displaystyle \sum _{k=1}^{n-1}{h}_k}\text{}\right)\right)\right\}\times }\\ \times \left\{\frac{\partial }{\partial M_m}\left({\alpha }_{l+1}{M}_{l+1}-{\alpha }_l{M}_l\right)-(1+i)\right.\times \\ \times \left.\left({\alpha }_{l+1}{M}_{l+1}-{\alpha }_l{M}_l\right)\frac{\partial }{\partial M_m}\left({\displaystyle \sum _{k=l+1}^n\frac{h_k}{{\delta }_k}}\right)\right\}.\end{array}$

(16)

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

 

$$\frac{\partial }{\partial M_m}\left({\alpha }_{l+1}{M}_{l+1}\text{}-{\alpha }_l{M}_l\right)=\left\{\begin{array}{c}{\alpha }_{l+1},l+1=m\\ -{\alpha }_l,l=m\end{array}\right..$$

(17)

$\frac{\partial }{\partial M_m}\left({\displaystyle \sum _{k=l+1}^n\frac{h_k}{{\delta }_k}}\right)=\left\{\begin{array}{c}-\frac{h_m}{2{\delta }_m{M}_m},l+1\le m\le n+1\\ \mathrm{0,}m<l+1\end{array}\right..$ (18)

Также,

$\begin{array}{c}\frac{\partial }{\partial {\kappa }_m}\left(\frac{\partial p_n(z)}{\partial z}\right)=\\ =\frac{i}{2{\delta }_n}{\epsilon }_0{e}^{i\omega t}{\displaystyle \sum _{l=0}^{n-1}\exp \left\{\left(1+i\right)\left(-{\displaystyle \sum _{k=l+1}^{n-1}\frac{h_k}{{\delta }_k}}-\frac{1}{{\delta }_n}\left(z-{\displaystyle \sum _{k=1}^{n-1}{h}_k}\right)\right)\right\}\times }\\ \times \left\{\left({\alpha }_{l+1}{M}_{l+1}-{\alpha }_l{M}_l\right)\frac{\partial }{\partial {\kappa }_m}\left({\displaystyle \sum _{k=l+1}^{n-1}\frac{h_k}{{\delta }_k}}\right)\right\},\end{array}$

(19)

где

$\frac{\partial }{\partial {\kappa }_m}\left({\displaystyle \sum _{k=l+1}^{n-1}\frac{h_k}{{\delta }_k}}\right)=\left\{\begin{array}{c}-\frac{h_m}{2{\delta }_m{\kappa }_m},l+1\le m\le n-1\\ \mathrm{0,}m<l+1\end{array}\right..$ (20)

Случай 2. Теперь обратимся к случаю, когда m = n, то есть, оценивается чувствительность вертикальной производной давления внутри некоторого слоя к его параметрам. В этом случае мы будем иметь:

$\begin{array}{c}\frac{\partial }{\partial M_m}\left(\frac{\partial p_n(z)}{\partial z}\right)=\\ =\frac{\left(1+i\right)}{2{\delta }_m}{\epsilon }_0{e}^{i\omega t}{\displaystyle \sum _{l=0}^{n-1}\exp \left\{\left(1+i\right)\left(-\text{}{\displaystyle \sum _{k=l+1}^{n-1}\frac{h_k}{{\delta }_k}}-\frac{1}{{\delta }_m}\left(z-\text{}{\displaystyle \sum _{k=1}^{n-1}{h}_k}\right)\right)\right\}\text{}\times }\\ \times \left\{\frac{\left({\alpha }_{l+1}{M}_{l+1}-{\alpha }_l{M}_l\right)}{2M_m}\left(1-\frac{(1+i)}{{\delta }_m}\left(z-{\displaystyle \sum _{k=1}^{n-1}{h}_k}\right)\right)\right.-\\ \left.-\frac{\partial }{\partial M_m}\left({\alpha }_{l+1}{M}_{l+1}-{\alpha }_l{M}_l\right)\right\}.\end{array}$

(21)

Здесь

$$\frac{\partial }{\partial M_m}\left({\alpha }_{l+1}{M}_{l+1}-{\alpha }_l{M}_l\right)=\left\{\begin{array}{c}{\alpha }_m,l+1=m\\ \mathrm{0,}l+1\ne m\end{array}\right..$$

(22)

Чувствительность к изменению коэффициента проницаемости определяется выражением:

$\begin{array}{c}\frac{\partial }{\partial {\kappa }_m}\left(\frac{\partial p_n(z)}{\partial z}\right)=\\ =\frac{\left(1+i\right)}{4{\delta }_m{\kappa }_m}{\epsilon }_0{e}^{i\omega t}{\displaystyle \sum _{l=0}^{n-1}\exp \left\{\left(1+i\right)\left(-\text{}{\displaystyle \sum _{k=l+1}^{n-1}\frac{h_k}{{\delta }_k}\text{}}-\frac{1}{{\delta }_n}\left(z-\text{}{\displaystyle \sum _{k=1}^{n-1}{h}_k}\right)\right)\right\}\text{}\times }\\ \times \left\{\left({\alpha }_{l+1}{M}_{l+1}-{\alpha }_l{M}_l\right)\left(1-(1+i)\frac{1}{{\delta }_m}\left(z-{\displaystyle \sum _{k=1}^{n-1}{h}_k}\right)\right)\right\}.\end{array}$

(23)

Полученные выражения позволяют проанализировать поведение чувствительности вертикальной компоненты электрического поля к изменению параметра α, модуля Био M и коэффициента проницаемости k в зависимости от глубины. Чувствительность к параметру αm слоя m равна нулю во всех вышележащих слоях, в то время как в слоях m ≤ l ≤ n она описывается одним и тем же выражением (с точностью до соответствующего сдвига по глубине) и имеет одинаковые максимальные значения на границах слоев (рис. 2а).

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

Сравнивая результаты расчета по формулам (19) и (23), рис. 2в, можно видеть, что внутри заданного слоя имеется чувствительность к коэффициенту проницаемости только самого этого слоя, в то время как чувствительность к проницаемости вышележащих слоев оказывается пренебрежимо мала.

 

Рис. 2. Графики чувствительности градиента порового давления к параметрам α (а), M (б) и k (в) второго слоя в трехслойной модели. Описание моделей приводится в табл. 1.

 

По соотношению графиков на рис. 2в можно видеть, что с ростом коэффициента k, чувствительность к его изменениям резко уменьшается. Одновременно происходит также и сокращение чувствительности к M (рис. 2б).

Резюмируя приведенные выше соображения, можно отметить, что градиент порового давления в n-слойной среде обладает чувствительностью к параметрам слоя m в диапазоне глубин, отвечающих слоям m ≤ l ≤ n.

3. ОСОБЕННОСТИ ОБРАТНОЙ ЗАДАЧИ ДЛЯ ГРАДИЕНТА ПОРОВОГО ДАВЛЕНИЯ В ГОРИЗОНТАЛЬНО-СЛОИСТОЙ СРЕДЕ. ЭКВИВАЛЕНТНОСТЬ И ПРЕДПОСЫЛКИ ПОВЫШЕНИЯ УСТОЙЧИВОСТИ

Рассматривая вопрос о решении обратной задачи, необходимо оговорить набор параметров, которые ищутся в процессе решения, а также тех параметров, которые считаются известными и задаются при решении прямой задачи на каждом шаге инверсии. Искомыми величинами являются коэффициент Био α, модуль Био M и проницаемость k, а к числу задаваемых параметров относятся вязкость порового флюида и объемная деформация ε₀, на которой мы остановимся несколько подробнее.

Общепринятое значение полусуточной приливной деформации 2 $‧$ 10^$–$8 [Латынина, Кармалеева, 1978]. Однако реальная (или измеряемая в скважине) объемная деформация ε₀ может отличаться от этого значения. В настоящей работе деформация считается заданной $(\epsilon =2\cdot {10}^{-8})$ лишь в целях упрощения и большей наглядности. На практике при решении обратной задачи деформацию можно считать неизвестной, объединив в общий параметр с коэффициентом Био, и определять произведение $\alpha {\epsilon }_0.$ Такой подход правомочен, поскольку обе названные величины входят в виде сомножителей только в правую часть уравнения для порового давления.

Поскольку величина коэффициента Био α представляет существенно меньший интерес по сравнению с модулем Био M и проницаемостью k, даже получение оценки произведения $\alpha {\epsilon }_0$ в совокупности с оценками M и k по отдельности можно считать достаточно важным результатом. Эти соображения никак не влияют на сделанные в работе выводы относительно эффективности применения предложенной методики для определения упругих характеристик и проницаемости формаций и ничтожно влияют на приведенную в работе схему решения обратной задачи.

Таким образом, даже допустив, что наличие скважины искажает деформации в спектре приливных колебаний [Segall et al., 2003], учет реальных (неизвестных) деформаций можно обеспечить, отказавшись от определения коэффициента Био и ограничившись оценкой более важных параметров: модуля Био и проницаемости. Тем не менее, в последующем изложении мы будем считать объемную деформацию известной и обсудим возможность определения всех трех характеристик: α, M и k.

Возвращаясь к выражениям (10) и (14), можно выделить несколько параметров, согласованные изменения которых сохраняют результирующие кривые p(z) и $\frac{\partial p(z)}{\partial z}.$ Очевидно, что любые изменения ${\alpha }_n,M_n,{\kappa }_n,$ сохраняющие величины ${\alpha }_n{M}_n$ и ${\delta }_n,$ не меняют решения в слое n и всех вышележащих слоях. Ясно, что широкий диапазон и логарифмический масштаб изменения параметров k и M могут приводить к их взаимокомпенсирующему влиянию и сохранению значения скин-толщины δ. Более того, из общих (гео)физических соображений следует, что такая ситуация весьма естественна, поскольку, например, увеличение упругого модуля горной породы согласуется с уменьшением коэффициента проницаемости.

Однако значительные изменения M приведут к таким изменениям произведения αM, которые не могут быть скомпенсированы соответствующими изменениями α, поскольку диапазон значений коэффициента Био весьма узок и составляет от 0 до 1 [Френкель, 1944].

Таким образом, область эквивалентных моделей, которым отвечают одинаковые решения

$p^{(1)}\left({\alpha }^{(1)},M^{(1)},{\kappa }^{(1)},{\eta }^{(1)}\right)=p^{(2)}\left({\alpha }^{(2)},M^{(2)},{\kappa }^{(2)},{\eta }^{(2)}\right),$

описывается следующими соотношениями:

 

$$\left\{\begin{array}{c}{\alpha }^{(1)}{M}^{(1)}={\alpha }^{(2)}{M}^{(2)}\\ \frac{{\kappa }^{(1)}{M}^{(1)}}{{\eta }^{(1)}}=\frac{{\kappa }^{(2)}{M}^{(2)}}{{\eta }^{(2)}}\end{array}\right..$$

(24)

Здесь индексы (1) и (2) отвечают двум эквивалентным моделям. Считая, что динамическая вязкость флюида не изменяется, система (24) сводится к условию:

 

$$\frac{{\alpha }^{(1)}}{{\kappa }^{(1)}}=\frac{{\alpha }^{(2)}}{{\kappa }^{(2)}}.$$

(25)

Другими словами, две модели со сколь угодно различающимися значениями M являются эквивалентными, если выполняется (25). На практике, однако, с учетом характерных диапазонов изменения α (0$–$1) и k (10$–$10$–$10$–$17), условие (25) описывает не слишком широкую область эквивалентности.

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

Как можно видеть, поровое давление и его градиент (а, соответственно, и вертикальное электрическое поле) в некотором слое зависят от параметров α, M и δ данного слоя, а также всех вышележащих слоев. Соответственно, значения, наблюденные в пределах этого слоя, будут содержать информацию о характеристиках всей толщи, за исключением нижележащей среды.

Это свойство существенно сужает область эквивалентности и снижает неустойчивость решения обратной задачи, позволяя восстановить параметры слоев по наблюденному (в скважине) вертикальному профилю $\frac{\partial p(z)}{\partial z}.$ Отметим, что набор мощностей h_k может быть легко получен на основе анализа $\frac{\partial p(z)}{\partial z}$ путем выделения точек (участков) разрыва (резкого изменения значений) последней. Соответственно, задавая значения h_k, можно избежать вышеуказанной неопределенности и разрешить модель по параметрам α, M и k в каждом слое.

4. АЛГОРИТМ 1-D ОБРАТНОЙ ЗАДАЧИ И РЕЗУЛЬТАТЫ ЕГО ОПРОБОВАНИЯ НА СИНТЕТИЧЕСКИХ НАБОРАХ ДАННЫХ

Для реконструкции этих параметров, на языке MATLAB с использованием функций библиотеки Optimization Toolbox была реализована процедура нелинейной инверсии. В качестве входных данных выступают значения вертикального градиента порового давления. Для решения прямой задачи используется выражение (14). Минимизируемый функционал представляет собой сумму квадратов абсолютных значений разностей $\frac{\partial p}{\partial z}\left(z_i\right)$ на сетке глубин z_i, шаг которой на порядок меньше минимальной скин-толщины δ_min:

 

$$\Phi ={\displaystyle \sum _{i=1}^N{\left|\frac{\partial p^{\text{obs}}}{\partial z}\left(z_i\right)-\frac{\partial p^{\mathrm{mod}}}{\partial z}\left(z_i,m\right)\right|}^2}.$$

(26)

Здесь: $\frac{\partial p^{\text{obs}}}{\partial z}$ $–$ наблюденные значения вертикального градиента порового давления (входные данные); $\frac{\partial p^{\mathrm{mod}}}{\partial z}$ $–$ теоретические значения, рассчитанные для модели, описываемой вектором m:

$m=\left\{{\alpha }_1\mathrm{,...,}{\alpha }_n,\lg M_1\mathrm{,...,}\lg M_n,\lg {\kappa }_1\mathrm{,...,}\lg {\kappa }_n\right\}.$ (27)

Алгоритм минимизации функционала (26) представляет собой модификацию квазиньютоновского метода [Powell, 1978; Coleman and Li, 1996]. При этом, подбор по параметрам k и M ведется в логарифмическом масштабе, по α $–$ в арифметическом.

Для оценки эффективности инверсии было проведено тестирование алгоритма на трех массивах синтетических распределений $\frac{\partial p}{\partial z}\left(z_i\right),$ рассчитанных для моделей из табл. 1 на равномерной сетке глубин в диапазоне от 0 до 500 м с шагом 5 м. Результаты решения обратной задачи представлены в табл. 2 в виде восстановленных значений искомых параметров, а также на рис. 3 в виде графиков входных данных и откликов, рассчитанных для результирующих моделей. Мощности слоев были заданы априорно. В качестве модели начального приближение выступала модель с одинаковыми значениями соответствующих параметров для всех трех слоев (см. второй столбец табл. 2). Учитывая отсутствие шумов во входных данных, в качестве критерия остановки процесса минимизации использовалась величина шага по параметрам исходной модели в некоторой норме.

Сопоставляя истинные значения параметров с полученными в процессе инверсии, можно видеть, что погрешность определения параметров k и M не превышает 10$–$15%, коэффициента α $–$ 20$–$25%. При этом, невязка между входными данными и откликом от результирующей модели составляет порядка десятых долей процента.

 

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

 

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

Параметр	Модель 1	Модель 2	Модель 3
h1, м	100	100	100
h2, м	150	150	150
α1	0.7	0.7	0.5
α2	0.5	0.5	0.7
α3	0.5	0.7	0.5
M1, Па	1010	1010	1010
M2, Па	109	109	109
M3, Па	1010	1010	1010
k1, м2	10–17	10–16	10–15
k2, м2	10–15	10–14	10–13
k3, м2	10–16	10–15	10–14
η1, Па ‧ с	10–3	10–3	10–3
η2, Па ‧ с	10–3	10–3	10–3
η3, Па ‧ с	10–3	10–3	10–3
δ1, м	1.17	3.71	11.7
δ2, м	3.71	11.7	37.1
δ3, м	3.71	11.7	37.1

 

Таблица 2. Результаты решения обратной задачи по набору синтетических данных

Параметр	Начальное приближение	Модель 1 (истинная)	Модель 1 (подбор)	Модель 2 (истинная)	Модель 2 (подбор)	Модель 3 (истинная)	Модель 3 (подбор)
α1	0.6	0.7	0.5962	0.7	0.5514	0.5	0.5120
α2	0.6	0.5	0.4466	0.5	0.4337	0.7	0.4532
α3	0.6	0.5	0.5321	0.7	0.5262	0.5	0.5074
M1	1011	1010	1010.0698	1010	1010.1036	1010	109.9897
M2	1011	109	109.0503	109	109.0618	109	109.1888
M3	1011	1010	109.9731	1010	1010.1240	1010	109.9936
k1	10–13	10–17	10–17.0697	10–16	10–16.1036	10–15	10–14.9897
k2	10–13	10–15	10–15.0502	10–14	10–14.0618	10–13	10–13.1888
k3	10–13	10–16	10–15.9731	10–15	10–15.1240	10–14	10–13.9936

 

5. ОБСУЖДЕНИЕ

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

О принципиальной возможности выделения приливных сигналов по измерениям в скважинах. На примере электрического каротажа методом собственных потенциалов (ПС) известно, что разности потенциалов, измеряемые между устьем скважины и заглубленным электродом, возникающие, в основном, вследствие электрохимических явлений в горных породах, составляют десятки и даже сотни милливольт [Serra, 1984]. Таким образом, выделение сигналов, порождаемых приливными деформациями и имеющих амплитуды порядка десятков и сотен микровольт на фоне более интенсивных (но слабо изменяющихся во времени) шумов электрохимического и иного происхождения, может оказаться нетривиальной задачей. Тем не менее, можно рассчитывать на ее успешное решение, принимая во внимание тот факт, что приливные колебания имеют известную (фиксированную) частоту и при достаточной длительности наблюдений данная узкополосная составляющая может быть выделена, например, путем накопления.

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

О существовании вертикального электрического поля электрокинетического происхождения на поверхности Земли. Как известно, электромагнитное поле, описываемое уравнениями Максвелла, в квазистационарном приближении не имеет вертикальной электрической составляющей на границе земля$–$воздух, если электропроводностью воздуха ? можно пренебречь:

Ez = 0 при z = 0, если σ(z) = 0

при z < 0, (28)

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

$E_z(0)=c\frac{\partial p}{\partial z}\left(0\right),$

где величина $\frac{\partial p}{\partial z}\left(0\right)$ отлична от нуля. Электрокинетический механизм возникновения поля представляет собой типичную стороннюю силу (сторонний возбудитель), чему соответствует дополнительное слагаемое в правой части уравнения Максвелла:

 

$$rotH=\sigma E+j_{ext},$$

(29)

где: $j_{ext}=\sigma (z)E_z^{ek}(z),$ а $E_z^{ek}(z)$ $–$ рассматриваемое нами электрическое поле электрокинетической природы. Идеализируя ситуацию, можно полагать, что поскольку в приземном слое (как и всюду ниже) имеется система вертикальных капиллярров, доходящих непосредственно до земной поверхности и давление поровой жидкости внутри каждого капилляра изменяется с глубиной, то между концами капилляра всегда имеется разность электрических потенциалов, пропорциональная разности давлений на них независимо от свойств вышележащей среды (воздуха). Другими словами, условие (28) не должно распространяться на поле $E_z^{ek}(z),$ существующее по причине наличия внешней силы (градиента порового давления), а не энергетических преобразований внутри чисто электродинамического процесса.

О переходе от электрического поля к полю порового давления. Основным соотношением теории электрокинетических явлений является формула Гельмгольца$–$Смолуховского:

 

$$E=c\cdot gradp,$$

(30)

где p $–$ поровое давление, а коэффициент

 

$$c=\frac{\epsilon \zeta }{4\pi \eta \sigma },$$

(31)

называют коэффициентом потокового потенциала. Здесь: ε $–$ диэлектрическая постоянная, ζ $–$ дзета-потенциал, η и σ $–$ вязкость и электропроводность флюида, соответственно. На практике величину потокового потенциала с в первом приближении можно считать постоянной, поскольку диапазон ее изменений (10^$–$6$–$10^$–$7 В/Па) значителен уже по сравнению с диапазоном изменений $\frac{\partial p}{\partial z}$ в пределах скин-слоя [Гершензон, Гохберг, 1994].

О зависимости изучаемых параметров от пористости. Несмотря на то, что и модуль Био M, и коэффициент проницаемости k зависят от коэффициента пористости $\varphi$, характер этой связи позволяет при решении обратной задачи рассматривать оба названных выше параметра как независимые. Зависимость проницаемости k от пористости $\varphi$ близка к экспоненциальной $(\log \kappa ~\varphi ),$ в то время как модуль Био в некотором приближении обратно пропорционален квадрату пористости $\left(M~\frac{1}{{\varphi }^2}\right).$ Таким образом, реакция проницаемости и модуля Био на изменение пористости оказывается существенно разной. Другими словами, при одном и том же изменении пористости, изменения проницаемости очень велики, а модуля Био $–$ несущественны, что позволяет извлекать информацию о пористости из их комбинации.

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

6. ЗАКЛЮЧЕНИЕ

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

Разработаны алгоритм и программа решения 1-D-обратной задачи, позволяющие восстанавливать значения $\alpha (z),M(z),\kappa (z)$ по известным распределениям $\frac{\partial p(z)}{\partial z}$ внутри горизонтально-слоистой среды. Оценка возможностей инверсии, проведенная на наборе синтетических данных, свидетельствует о достаточно высокой точности реконструкции упругих параметров и коэффициентов проницаемости слоев. Необходимо отметить, что на практике, при выполнении инверсии по электрическому полю, величину потокового потенциала с в первом приближении можно считать постоянной (или кусочно-постоянной), вводя в качестве дополнительного параметра в алгоритм обратной задачи.

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

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

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

Настоящая работа является развитием исследований, начатых при непосредственном участии И.А. Гарагаша и Н.И. Колосницына, которым авторы выражают свою искреннюю признательность за содержательные консультации и поддержку. Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект 17-05-00511_а) и Госзадания ИФЗ РАН.

About the authors

D. A. Alekseev

Institute of the Earth Physics of the Russian Academy of Sciences; P.P. Shirshov Institute of Oceanology of Russian Academy of Sciences; Moscow Institute of Physics and Technology
 

Author for correspondence.
Email: alexeevgeo@gmail.com
Russian Federation, Bolshaya Gruzinskaya str., 10-1, Moscow 123242, Russia; 36, Nakhimovskii prospect, Moscow, 117997; 9, Institutskij, Dolgoprudny, Moscow region, 141701

M. B. Gokhberg

Institute of the Earth Physics of the Russian Academy of Sciences

Email: gmb@ifz.ru
Russian Federation, Bolshaya Gruzinskaya str., 10-1, Moscow 123242, Russia

References

Supplementary files