Физика Земли, 2022, № 3, стр. 92-107

ВВЕДЕНИЕ

Решение прямых и обратных трехмерных задач геоэлектрики бианизотропных сред представляет актуальную проблему современной электроразведки, поскольку такие среды позволяют максимально полно описать линейные электромагнитные свойства горной породы [Туров, 1983; Александров, 2000]. Сложность решения таких задач связана с большим количеством электромагнитных параметров (максимально 36), которыми характеризуются линейные среды и сложностью решения уравнений относительно компонент электромагнитного поля для таких сред. Одним из способов решения прямых задач для бианизотропных сред является способ сведения системы уравнений Максвелла к системе обыкновенных дифференциальных уравнений первого порядка. Однако такой подход эффективен только для одномерных сред. Для трехмерных сред наиболее перспективным способом является способ, основанный на решении объемных интегральных уравнений. Вследствие этого, прежде получим объемные интегральные уравнения относительно компонент электромагнитного поля для случая трехмерной бианизотропной неоднородности, находящейся в однородном бианизотропном пространстве.

ВЫВОД ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ

Для решения прямой трехмерной задачи геоэлектрики бианизотропных сред воспользуемся методом интегральных уравнений на основе леммы Лоренца [Марков и др., 1979].

Рассмотри две системы уравнений Максвелла:

где: $i = \sqrt { - 1} $; ${\mathbf{H}}$, ${\mathbf{E}}$ – векторы напряженностей соответственно магнитного и электрического полей, источниками которых являются компоненты стороннего электрического тока ${{{\mathbf{J}}}^{{ext}}}$ и сторонней магнитной индукции ${{{\mathbf{B}}}^{{ext}}}$ (в уравнениях (1) и (2)); ${\mathbf{G}}_{x}^{e}$, ${\mathbf{Q}}_{x}^{e}$ – вспомогательные векторы, аналогичные напряженностям соответственно магнитного и электрического полей от источников, сосредоточенных в точках, направленными по соответствующим осям координат возбудителями электромагнитного поля в виде пространственных дельта – функций Дирака (уравнение (3)); ${\mathbf{i}}$ – орт оси $x$ декартовой системы координат; $\omega $ – частота; $\sigma $, $\mu $, $\alpha $, $\beta $ – электромагнитные параметры неоднородной среды; ${{\sigma }_{0}}$, ${{\mu }_{0}}$, ${{\alpha }_{0}}$, ${{\beta }_{0}}$ – электромагнитные параметры вмещающей среды.

Магнитная проницаемость $\mu $ и удельная электропроводность $\sigma $ – известные параметры в геоэлектрике, $\alpha $, $\beta $ – бианизотропные параметры, физическая размерность которых, соответственно, равны (${1 \mathord{\left/ {\vphantom {1 {\text{м}}}} \right. \kern-0em} {\text{м}}}$) и (${1 \mathord{\left/ {\vphantom {1 {\text{м}}}} \right. \kern-0em} {\text{м}}} \cdot {\text{c}}$). Все эти параметры описывают наиболее общие линейные электромагнитные свойства среды и являются матрицами размерностью 3 × 3 и, следовательно, общее количество электромагнитных параметров равно 36. Как показано в работе [Александров, 2000], бианизотропные параметры $\alpha $, $\beta $ связаны с извилистостью токопроводящих путей. Они являются комплексными функциями частоты и равны нулю при $\omega = 0$, т.е. при постоянном токе они исчезают.

Отметим, что исследованиями бианизотропных сред занимался один из основоположников теории геоэлектрики В.Р. Бурсиан [Bursian, Timorew, 1926; Бурсиан, 1972], получивший фундаментальные научные результаты, имеющие огромное значение и в настоящее время. Рассмотрение уравнений Максвелла в наиболее общей линейной постановке связано с объектом геофизических исследований, для которого заранее нельзя определить конкретный вид материальных уравнений.

Умножим уравнение (1) на транспонированный вектор электрического поля ${\mathbf{Q}}{{_{x}^{e}}^{T}}$, уравнение (4) на транспонированный вектор магнитного поля ${{{\mathbf{H}}}^{T}}$, и результаты умножения вычтем:

Далее, умножим уравнение (3) на ${{{\mathbf{E}}}^{T}}$, уравнение (2) на транспонированный вектор электрического поля ${\mathbf{G}}_{x}^{T}$ и результаты умножения вычтем:

Вычитая результаты этих операций, получим:

Для остальных компонент электрического поля (меняя ориентацию источников в уравнении (3) на $y$-компоненту с ортом ${\mathbf{j}}$ и $z$-компоненту с ортом ${\mathbf{k}}$), действуя аналогично, получим:

Отсюда, после интегрирования по всему пространству, получим:

Аналогичные интегральные уравнения получим, введя точечный источник сторонней магнитной индукции:

Тогда:

Вводя вектор ${\mathbf{X}} = \left( {\begin{array}{*{20}{c}} {\mathbf{E}} \\ {\mathbf{H}} \end{array}} \right)$, получим векторные объемные интегральные уравнения в случае бианизотропных сред:

где ${{{\mathbf{X}}}^{f}} = \int_{{{V}_{\infty }}} {\left( {\begin{array}{*{20}{c}} {{{{\hat {Q}}}_{e}}}&{{{{\hat {G}}}_{e}}} \\ {{{{\hat {Q}}}_{h}}}&{{{{\hat {G}}}_{h}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{{{\mathbf{J}}}^{{ext}}}} \\ {{{{\mathbf{B}}}^{{ext}}}} \end{array}} \right)dv} $.

Здесь интегрирование производится по области источников и носителей неоднородностей во вмещающей среде с параметрами ${{\sigma }_{0}}$, ${{\mu }_{0}}$, ${{\alpha }_{0}}$, ${{\beta }_{0}}$.

РЕШЕНИЕ ПРЯМОЙ ЗАДАЧИ НА ОСНОВЕ ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ

Для практических вычислений в качестве вмещающей среды использовалось однородное изотропное пространство. Найдем выражения для тензорной функции Грина электрического ${{\hat {Q}}_{e}} = \left( {\begin{array}{*{20}{c}} {{\mathbf{Q}}{{{_{x}^{e}}}^{T}}} \\ {{\mathbf{Q}}{{{_{y}^{e}}}^{T}}} \\ {{\mathbf{Q}}{{{_{z}^{e}}}^{T}}} \end{array}} \right)$ и магнитного ${{\hat {G}}_{e}} = \left( {\begin{array}{*{20}{c}} {{\mathbf{G}}_{x}^{{eT}}} \\ {{\mathbf{G}}_{y}^{{eT}}} \\ {{\mathbf{G}}_{z}^{{eT}}} \end{array}} \right)$ типов. Ротируя систему уравнений (3) и (4) и учитывая ${\text{div}}{\mathbf{G}}_{x}^{e} = 0$ и ${\text{div}}{\mathbf{Q}}_{x}^{e} = - \frac{1}{{{{\sigma }_{0}}}}{\text{div}}\delta {\mathbf{i}}$, получим:

Тогда:

где $G = - \frac{1}{{4\pi R}}{{e}^{{ - \sqrt {i\omega \mu {{\sigma }_{0}}} R}}}$, $\operatorname{Re} (\sqrt {i\omega \mu {{\sigma }_{0}}} ) > 0$, $R = \sqrt {{{{(x - x{\kern 1pt} ')}}^{2}} + {{{(y - y{\kern 1pt} ')}}^{2}} + {{{(z - z{\kern 1pt} ')}}^{2}}} $ есть фундаментальное обобщенное решение уравнения Гельмгольца:

Учитывая, что ${{\hat {Q}}_{h}} = {{\hat {G}}_{e}}$, ${{\hat {G}}_{h}} = - \frac{{{{\sigma }_{0}}}}{{i\omega {{\mu }_{0}}}}{{\hat {Q}}_{e}}$, получим:

или

Последнее уравнение можно привести к виду:

где и

$Q = \left( {\begin{array}{*{20}{c}} {{{{\hat {Q}}}_{e}}}&{{{{\hat {G}}}_{e}}} \\ {{{{\hat {G}}}_{e}}}&{ - \frac{{{{\sigma }_{0}}}}{{i\omega {{\mu }_{0}}}}{{{\hat {Q}}}_{e}}} \end{array}} \right);$ $G = \left( {\begin{array}{*{20}{c}} {{{{\hat {Q}}}_{e}}}&{{{{\hat {G}}}_{e}}} \\ {{{{\hat {G}}}_{e}}}&{ - \frac{{{{\sigma }_{0}}}}{{i\omega {{\mu }_{0}}}}{{{\hat {Q}}}_{e}}} \end{array}} \right)$ – тензорные функции Грина электродинамической задачи. Теория функций Грина хорошо разработана (см., например, [Дмитриев, 2012]), поэтому в настоящей работе явный вид этих величин выписываться не будет.

Полагая, что в области элементарного объема $dV = \Delta x\Delta y\Delta z$ поле и электромагнитные параметры среды постоянны, переход к системе линейных алгебраических уравнений будет связан с вычислением интегралов по элементарному объему от тензорных функций Грина. Для тензорной функции Грина электрического типа получим:

где:

Для тензорной функции Грина магнитного типа получим:

РЕШЕНИЕ ОБРАТНОЙ ЗАДАЧИ ГЕОЭЛЕКТРИКИ НА ОСНОВЕ ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ

Решения прямых задач для большинства физических полей, используемых при изучении геологической среды, могут быть сведены к решению интегральных уравнений вида (5), где ${\mathbf{X}}$ – вектор напряженностей электромагнитного поля; ${{{\mathbf{X}}}^{f}}(x,y,z,\omega )$ – вектор первичного поля – решение прямой задачи для среды, для которой имеется функция Грина $G(x{\kern 1pt} ',y{\kern 1pt} ',z{\kern 1pt} ',x,y,z,\omega )$; $\Delta S(x,y,z) = $ $S(x,y,z) - {{S}^{0}}$ – избыточные значения физической величины неоднородностей, которые являются искомыми, где ${{S}^{0}}$ – электромагнитные параметры вмещающей среды, $S(x,y,z)$ – электромагнитные параметры неоднородностей.

В дискретном виде, разбивая интеграл на сумму в выражении (5), получим:

где индексы означают: $p$ – точку наблюдения; $s$ – точку источника; $\nu $ – точку в неоднородности; $\Delta S$ – квадратная матрица электромагнитных параметров, имеющая квазидиагональный (принципиально для дальнейших вычислений) вид: где: подматрицы $\Delta {{S}_{j}}$ есть тензоры избыточных электромагнитных параметров (матрица размерностью 6 × 6 элементов – в общем случае, 36 параметров) $j$-го элемента объема всей неоднородности $V$, состоящей из $N$ элементов; ${\mathbf{X}}_{p}^{s}$ – первичное поле от источника в точке наблюдения; ${{{\mathbf{X}}}_{v}}$ – поле в неоднородности; ${{{\mathbf{X}}}^{p}}$ – поле в точке наблюдения (как разность потенциалов); $G_{v}^{p}$ – передаточная матрица от неоднородности в точку наблюдения, размерности $K \times 6N$; $K$ – количество точек наблюдения.

Найдем поле в неоднородностях ${{{\mathbf{X}}}_{v}} = $ $ = G_{v}^{v}\Delta S{{{\mathbf{X}}}_{v}} + {\mathbf{X}}_{v}^{s}$, тогда ${{{\mathbf{X}}}_{v}} = {{([1] - G_{v}^{v}\Delta S)}^{{ - 1}}}{\mathbf{X}}_{v}^{s}$.

Отсюда ${{{\mathbf{X}}}^{p}} = G_{v}^{p}\Delta S{{{\mathbf{X}}}_{v}} + {\mathbf{X}}_{p}^{s} = $ $G_{v}^{p}\Delta S([1] - $ $ - \,\,G_{v}^{v}\Delta S{{)}^{{ - 1}}}{\mathbf{X}}_{v}^{s} + {\mathbf{X}}_{p}^{s}$.

Здесь ${\mathbf{X}}_{v}^{s}$ – поле от источника в неоднородность, вычисляется через функцию Грина вмещающей среды, ${\mathbf{X}}_{p}^{s}$ – первичное поле в точке наблюдения, вычисляется через функцию Грина вмещающей среды.

Или для конкретного источника

где $k$ – номер источника и воспользовались правилом $(A{{B}^{{ - 1}}}) = {{(B{{A}^{{ - 1}}})}^{{ - 1}}}$, где$A$ и $B$ – некоторые неособенные квадратные матрицы.

Используя дополнительные источники, введем составные матрицы:

где $M$ – количество источников.

Перейдем от векторной системы уравнений к матричной:

Для переопределенной и нормально определенной системы из цепочки уравнений:

получим:

Перепишем последнее уравнение в виде $\Delta {{S}^{{ - 1}}}A = B,$ где $A = {{(G{{_{v}^{p}}^{T}}G_{v}^{p})}^{{ - 1}}}G{{_{v}^{p}}^{T}}(X - Z) = \{ {{a}_{{ij}}}\} $, $B = Y + $ $ + \,\,G_{v}^{v}{{(G{{_{v}^{p}}^{T}}G_{v}^{p})}^{{ - 1}}}G{{_{v}^{p}}^{T}}(X - Z) = \{ {{b}_{{ij}}}\} $.

Последнее уравнение с учетом представления (6) разбивается на отдельные уравнения по следующе схеме:

Отсюда нахождение искомых параметров разбивается на решение отдельных систем матричных уравнений

Подматрицы ${{a}_{{ij}}}$ имеют размеры 6 × 1 (вектор-столбец). Для разрешимости каждого из этих уравнений необходимо минимум 6 источников. В случае 6 источников система уравнений является нормально определенной, в случае большего количества источников система уравнений становиться переопределенной, в случае менее 6 источников система уравнений становится недоопределенной и в дальнейшем последний случай рассматриваться не будет.

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

ПРИМЕР РЕШЕНИЯ ПРЯМОЙ И ОБРАТНОЙ ЗАДАЧИ НА ОСНОВЕ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ

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

С использованием изложенного подхода был разработан алгоритм и создано программное обеспечение математического моделирования электромагнитных полей в линейных трехмерных неоднородных средах, при этом, в целях упрощения изложения, функция Грина вмещающей среды взята для однородного пространства. Программа позволяет вычислять все компоненты электромагнитного поля ${{E}_{x}}$, ${{E}_{y}}$, ${{E}_{z}}$, ${{H}_{x}}$, ${{H}_{y}}$, ${{H}_{z}}$ от всех возможных видов источников – ${{J}_{x}}^{{ext}}$, ${{J}_{y}}^{{ext}}$, ${{J}_{z}}^{{ext}}$, ${{B}_{x}}^{{ext}}$, ${{B}_{y}}^{{ext}}$, ${{B}_{z}}^{{ext}}$. Количество сторонних источников электрического типа (${{J}_{x}}^{{ext}}$, ${{J}_{y}}^{{ext}}$, ${{J}_{z}}^{{ext}}$), и сторонних источников магнитного типа (${{B}_{x}}^{{ext}}$, ${{B}_{y}}^{{ext}}$, ${{B}_{z}}^{{ext}}$), сосредоточенных в одной точке равно шести ($k = 6$). Этого достаточно, чтобы сформировать составные матрицы в выражении (7) и разрешить системы линейных алгебраических уравнений (9) относительно искомых электромагнитных параметров неоднородностей.

Рассмотрим изотропную вмещающую среду с параметрами воздуха: ${{\sigma }_{0}} = {{10}^{{ - 14}}} + i\omega \varepsilon ,$ $\varepsilon = \frac{1}{{36\pi }}{{10}^{{ - 9}}}\,\,{{\text{Ф}} \mathord{\left/ {\vphantom {{\text{Ф}} {\text{м}}}} \right. \kern-0em} {\text{м}}}$, ${{\mu }_{0}} = 4\pi {{10}^{{ - 7}}}\,\,{{{\text{Гн}}} \mathord{\left/ {\vphantom {{{\text{Гн}}} {\text{м}}}} \right. \kern-0em} {\text{м}}}$, ${{\alpha }_{0}} = {{\beta }_{0}} = [0]$, где $[0]$ – нулевая матрица размерностью 3 × 3.

Рассмотрим объемный источник, частотой $\omega = 2\pi {{10}^{6}}\,{\text{Гц}}$, в виде куба со сторонами $\Delta x = 0.01\,{\text{м}}$, $\Delta y = 0.01\,{\text{м}}$, $\Delta z = 0.01\,{\text{м}}$, параллельными осям координат, центр которого находится в точке с координатами (0,0, ${{Z}_{s}} = 0.1\,{\text{м}}$). На рис. 1 его местоположение отмечено окружностью черного цвета на соответствующей глубине. Местоположение приемников в плоскости $Z = 0$ обозначены красными точками, неоднородность – звездочкой синего цвета, которая имеет размер $0.1\,{\text{м}} \times 0.1\,{\text{м}} \times 0.1\,{\text{м}}$ и находится на глубине $2\,{\text{м}}$ (рис.1).

Пример 1. Решение прямой задачи для изотропной неоднородности.

Неоднородность – проводящий объект с параметрами

где $[1]$ – единичная матрица размерностью 3 × 3.

На рис. 2 приведены только действительная и мнимая часть $x$ – компоненты напряженности электрического поля ${{E}_{x}}$ от источника $J_{x}^{{ext}}$.

Пример 2. Решение прямой задачи для бианизотропной неоднородности.

Неоднородность – проводящий объект с параметрами

На рис. 3 приведены действительные и мнимые части всех компонент электромагнитного поля от всех источников.

Пример 3. Решение обратной задачи на основе примера 2.

Ниже приведены результаты решения обратной задачи на основе изложенного выше алгоритма. Индекс “исх” означает исходный тензор, “найд” – найденные значения тензора на основе алгоритма решения обратной задачи. В таблицах приведено покомпонентное сравнение величин.

Тензор удельной электропроводности:

Относительная погрешность ${{\left\| {\sigma _{{}}^{{{\text{исх}}}} - \sigma _{{}}^{{{\text{найд}}}}} \right\|} \mathord{\left/ {\vphantom {{\left\| {\sigma _{{}}^{{{\text{исх}}}} - \sigma _{{}}^{{{\text{найд}}}}} \right\|} {\left\| {\sigma _{{}}^{{{\text{исх}}}}} \right\|}}} \right. \kern-0em} {\left\| {\sigma _{{}}^{{{\text{исх}}}}} \right\|}} \times 100\% = 0.0263\% .$

Тензор магнитной проницаемости:

Относительная погрешность ${{\left\| {\mu _{{}}^{{{\text{исх}}}} - \mu _{{}}^{{{\text{найд}}}}} \right\|} \mathord{\left/ {\vphantom {{\left\| {\mu _{{}}^{{{\text{исх}}}} - \mu _{{}}^{{{\text{найд}}}}} \right\|} {\left\| {\mu _{{}}^{{{\text{исх}}}}} \right\|}}} \right. \kern-0em} {\left\| {\mu _{{}}^{{{\text{исх}}}}} \right\|}} \times 100\% = 0.0378\% .$

Бианизотропный параметр $\alpha $:

Относительная погрешность решения обратной задачи относительно бианизотропного параметра $\alpha $ ${{\left\| {\alpha _{{}}^{{{\text{исх}}}} - \alpha _{{}}^{{{\text{найд}}}}} \right\|} \mathord{\left/ {\vphantom {{\left\| {\alpha _{{}}^{{{\text{исх}}}} - \alpha _{{}}^{{{\text{найд}}}}} \right\|} {\left\| {\alpha _{{}}^{{{\text{исх}}}}} \right\|}}} \right. \kern-0em} {\left\| {\alpha _{{}}^{{{\text{исх}}}}} \right\|}} \times 100\% = 0.0301\% .$

Бианизотропный параметр $\beta $:

Относительная погрешность решения обратной задачи относительно бианизотропного параметра $\beta $ составляет ${{\left\| {\beta _{{}}^{{{\text{исх}}}} - \beta _{{}}^{{{\text{найд}}}}} \right\|} \mathord{\left/ {\vphantom {{\left\| {\beta _{{}}^{{{\text{исх}}}} - \beta _{{}}^{{{\text{найд}}}}} \right\|} {\left\| {\beta _{{}}^{{{\text{исх}}}}} \right\|}}} \right. \kern-0em} {\left\| {\beta _{{}}^{{{\text{исх}}}}} \right\|}} \times 100\% = 0.0215\% .$

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

ВЫВОДЫ

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

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

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