Физика Земли, 2022, № 2, стр. 136-143

ВВЕДЕНИЕ

Обратную задачу геофизики можно разделить на две последовательные подзадачи: 1) определение местоположения и геометрии неоднородности как аномалиеобразующего объекта; и 2) определение физических свойств горных пород, составляющих этот объект.

Первая подзадача связана с истокообразной аппроксимацией, результатом которой является обнаружение аномалиеобразующих тел и определение их “геометрии” [Александров, Монахов, 2014]. Этот вопрос в настоящей работе рассматриваться не будет. Вторая задача заключается в насыщении аномалиеобразующих объектов физическими параметрами. Этому вопросу посвящена настоящая статья.

ОСНОВНАЯ ИДЕЯ М.В. КЛИБАНОВА РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ ГЕОРАДАРНОГО МЕТОДА, ИЗЛОЖЕННАЯ В РАБОТЕ [BEILINA, KLIBANOV, 2012]

Пусть в среде возбуждается, к примеру, электромагнитное поле, удовлетворяющее дифференциальному уравнению в частных производных второго порядка, получаемому на основе уравнений Максвелла [Максвелл, 1989]. Задача сводится к уравнению Гельмгольца вида ${{\nabla }^{2}}\varphi - {{(i\omega )}^{n}}\mu \sigma \varphi = 0$, где φ – физическое поле, $\omega $ – частота, σ – удельная электропроводность, μ – магнитная проницаемость вакуума, $i = \sqrt { - 1} $ – мнимая единица.

Выразим неизвестную искомую функцию удельной электропроводности среды, преобразуя исходное уравнение к виду:

Отметим здесь, что функция электропроводности материальной среды $\sigma (x,y,z)$ не зависит от частоты возбуждаемого поля. После дифференцирования обеих частей равенства (1) по частоте, получим нелинейное дифференциальное уравнение, не содержащее электропроводность $\sigma (x,y,z)$:

Решая краевую задачу для нелинейного дифференциального уравнения (2), найдем поле $\varphi $. Подставляя его в уравнение (1), получим в пространстве искомое распределение электропроводности $\sigma (x,y,z)$. Удельная электропроводность рассматривается здесь как скалярная функция пространственных координат.

Покажем, что данный подход может быть распространен и на более сложные модели.

Из электродинамики сплошных сред известно, что для решения уравнений Максвелла, связывающих электрическое ${\mathbf{E}}$ и магнитное ${\mathbf{H}}$ поля через систему векторных дифференциальных уравнений первого порядка, которые в частотной области имеют вид:

где ${\mathbf{J}}$ – плотность электрического тока; ${\mathbf{B}}$ – индукция магнитного поля; ${{{\mathbf{I}}}^{{ex}}}$ – плотность контролируемого стороннего электрического тока, необходимо установить связи ${\mathbf{J}} = {\mathbf{J}}({\mathbf{E}},{\mathbf{H}})$ и ${\mathbf{B}} = {\mathbf{B}}({\mathbf{E}},{\mathbf{H}})$, т.е. определить материальные соотношения. Определив эти соотношения, уравнения Максвелла приобретают замкнутую форму, что дает возможность решения прямых и обратных задач.

Параметры среды являются коэффициентами, которые входят в материальные уравнения поля. В наиболее общем виде они могут быть определены в любой точке пространства разложением в ряд Тейлора по малым величинам напряженностей полей ${\mathbf{E}}$ и ${\mathbf{H}}$, что и имеет место в геоэлектрике [Туров, 1983]. Так, для плотности электрического тока получим:

где ${{\nabla }^{e}} = {\mathbf{i}}\frac{\partial }{{\partial {{E}_{x}}}} + {\mathbf{j}}\frac{\partial }{{\partial {{E}_{y}}}} + {\mathbf{k}}\frac{\partial }{{\partial {{E}_{z}}}}$, ${{\nabla }^{h}} = {\mathbf{i}}\frac{\partial }{{\partial {{H}_{x}}}} + {\mathbf{j}}\frac{\partial }{{\partial {{H}_{y}}}} + $ $ + \,\,{\mathbf{k}}\frac{\partial }{{\partial {{H}_{z}}}}$, причем производные находятся при ${\mathbf{E}},{\mathbf{H}} = 0$.

Матрица коэффициентов $\hat {\sigma } = \left( {\begin{array}{*{20}{c}} {\frac{{\partial {{J}_{x}}}}{{\partial {{E}_{x}}}}}&{\frac{{\partial {{J}_{x}}}}{{\partial {{E}_{y}}}}}&{\frac{{\partial {{J}_{x}}}}{{\partial {{E}_{z}}}}} \\ {\frac{{\partial {{J}_{y}}}}{{\partial {{E}_{x}}}}}&{\frac{{\partial {{J}_{y}}}}{{\partial {{E}_{y}}}}}&{\frac{{\partial {{J}_{y}}}}{{\partial {{E}_{z}}}}} \\ {\frac{{\partial {{J}_{z}}}}{{\partial {{E}_{x}}}}}&{\frac{{\partial {{J}_{z}}}}{{\partial {{E}_{y}}}}}&{\frac{{\partial {{J}_{z}}}}{{\partial {{E}_{z}}}}} \end{array}} \right) = $ $ = \left( {\begin{array}{*{20}{c}} {{{\sigma }_{{11}}}}&{{{\sigma }_{{12}}}}&{{{\sigma }_{{13}}}} \\ {{{\sigma }_{{21}}}}&{{{\sigma }_{{22}}}}&{{{\sigma }_{{23}}}} \\ {{{\sigma }_{{31}}}}&{{{\sigma }_{{32}}}}&{{{\sigma }_{{33}}}} \end{array}} \right)$ определяет тензор электропроводности в линейном законе Ома. Элементы этого тензора, как и другие электромагнитные параметры, могут зависеть от частоты $\omega $. Следующие слагаемые, связанные с оператором ${{\nabla }^{e}}$, определяют нелинейные электрические свойства среды. В силу малости напряженностей электромагнитного поля в геологической среде ими можно пренебречь, и в дальнейшем рассматривать только линейные материальные соотношения.

Матрица коэффициентов $\hat {\xi } = \left( {\begin{array}{*{20}{c}} {\frac{{\partial {{J}_{x}}}}{{\partial {{H}_{x}}}}}&{\frac{{\partial {{J}_{x}}}}{{\partial {{H}_{y}}}}}&{\frac{{\partial {{J}_{x}}}}{{\partial {{H}_{z}}}}} \\ {\frac{{\partial {{J}_{y}}}}{{\partial {{H}_{x}}}}}&{\frac{{\partial {{J}_{y}}}}{{\partial {{H}_{y}}}}}&{\frac{{\partial {{J}_{y}}}}{{\partial {{H}_{z}}}}} \\ {\frac{{\partial {{J}_{z}}}}{{\partial {{H}_{x}}}}}&{\frac{{\partial {{J}_{z}}}}{{\partial {{H}_{y}}}}}&{\frac{{\partial {{J}_{z}}}}{{\partial {{H}_{z}}}}} \end{array}} \right) = $ $ = \left( {\begin{array}{*{20}{c}} {{{\xi }_{{11}}}}&{{{\xi }_{{12}}}}&{{{\xi }_{{13}}}} \\ {{{\xi }_{{21}}}}&{{{\xi }_{{22}}}}&{{{\xi }_{{23}}}} \\ {{{\xi }_{{31}}}}&{{{\xi }_{{32}}}}&{{{\xi }_{{33}}}} \end{array}} \right)$ определяет тензор бианизотропных параметров.

Аналогично получим формулы и для магнитной индукции:

Матрица коэффициентов $\hat {\mu } = \left( {\begin{array}{*{20}{c}} {\frac{{\partial {{B}_{x}}}}{{\partial {{H}_{x}}}}}&{\frac{{\partial {{B}_{x}}}}{{\partial {{H}_{y}}}}}&{\frac{{\partial {{B}_{x}}}}{{\partial {{H}_{z}}}}} \\ {\frac{{\partial {{B}_{y}}}}{{\partial {{H}_{x}}}}}&{\frac{{\partial {{B}_{y}}}}{{\partial {{H}_{y}}}}}&{\frac{{\partial {{B}_{y}}}}{{\partial {{H}_{z}}}}} \\ {\frac{{\partial {{B}_{z}}}}{{\partial {{H}_{x}}}}}&{\frac{{\partial {{B}_{z}}}}{{\partial {{H}_{y}}}}}&{\frac{{\partial {{B}_{z}}}}{{\partial {{H}_{z}}}}} \end{array}} \right) = $ $ = \left( {\begin{array}{*{20}{c}} {{{\mu }_{{11}}}}&{{{\mu }_{{12}}}}&{{{\mu }_{{13}}}} \\ {{{\mu }_{{21}}}}&{{{\mu }_{{22}}}}&{{{\mu }_{{23}}}} \\ {{{\mu }_{{31}}}}&{{{\mu }_{{32}}}}&{{{\mu }_{{33}}}} \end{array}} \right)$ определяет тензор магнитной проницаемости. Элементы этого тензора также могут зависеть от частоты $\omega $.

Матрица коэффициентов $\hat {\zeta } = \left( {\begin{array}{*{20}{c}} {\frac{{\partial {{B}_{x}}}}{{\partial {{E}_{x}}}}}&{\frac{{\partial {{B}_{x}}}}{{\partial {{E}_{y}}}}}&{\frac{{\partial {{B}_{x}}}}{{\partial {{E}_{z}}}}} \\ {\frac{{\partial {{B}_{y}}}}{{\partial {{E}_{x}}}}}&{\frac{{\partial {{B}_{y}}}}{{\partial {{E}_{y}}}}}&{\frac{{\partial {{B}_{y}}}}{{\partial {{E}_{z}}}}} \\ {\frac{{\partial {{B}_{z}}}}{{\partial {{E}_{x}}}}}&{\frac{{\partial {{B}_{z}}}}{{\partial {{E}_{y}}}}}&{\frac{{\partial {{B}_{z}}}}{{\partial {{E}_{z}}}}} \end{array}} \right) = $ $ = \left( {\begin{array}{*{20}{c}} {{{\zeta }_{{11}}}}&{{{\zeta }_{{12}}}}&{{{\zeta }_{{13}}}} \\ {{{\zeta }_{{21}}}}&{{{\zeta }_{{22}}}}&{{{\zeta }_{{23}}}} \\ {{{\zeta }_{{31}}}}&{{{\zeta }_{{32}}}}&{{{\zeta }_{{33}}}} \end{array}} \right)$ является параметром бианизотропии и в совокупности с $\hat {\xi }$ определяет бианизотропные свойства среды [Bursian, Timorew, 1926; Третьяков, 1994].

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

Слагаемые ${{{\mathbf{J}}}^{{ex}}} = {\mathbf{J}}({\mathbf{0}},{\mathbf{0}})$ и ${{{\mathbf{B}}}^{{ex}}} = {\mathbf{B}}({\mathbf{0}},{\mathbf{0}})$ определяют активные свойства среды. При равенстве нулю сторонних токов они, в общем случае, не пропадают и могут определять неконтролируемые источники электромагнитного поля, распределенные в геологической среде, что является основой пассивного электромагнитного мониторинга современных геодинамических процессов [Александров и др., 2018].

Таким образом, наиболее общие линейные свойства геоэлектрической среды в частотной области отражаются в следующих материальных уравнениях, связывающих плотность электрического тока ${\mathbf{J}}$ и магнитную индукцию ${\mathbf{B}}$ с напряженностями электрического ${\mathbf{E}}$ и магнитного ${\mathbf{H}}$ полей:

где $\hat {\sigma }$, $\hat {\xi }$, $\hat {\mu }$, $\hat {\zeta }$ – являются матрицами коэффициентов размерности 3 × 3, зависящих от частоты $\omega $, причем $\hat {\sigma }$, $\hat {\mu }$ – соответственно, матрицы удельной электропроводности и магнитной проницаемости среды, хорошо известные в электродинамике. Величины $\hat {\zeta }$, $\hat {\xi }$ являются новыми в практике геофизики. Они связаны со сложной геометрией проводящих капилляров [Александров, 2000].

Отсюда следует, что общее количество электромагнитных параметров, описывающих наиболее общие линейные свойства геоэлектрической среды, равно 36. Отметим, что в законе Гука теории упругости, упругих параметров – 81.

Несмотря на большое количество физических параметров, подход М.В. Клибанова может быть использован для решения обратных задач и для полученных выше материальных уравнений.

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

Решения прямых задач для большинства физических полей, используемых при изучении геологической среды, могут быть сведены к решению интегральных уравнений вида [Морс, Фешбах, 1960]:

где ${\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)$ – избыточные значения физической величины неоднородностей, которые являются искомыми.

Для определенности будем рассматривать постоянный электрической ток как наиболее простую модель потенциального физического поля. Тогда $\Delta S(x,y,z) = S(x,y,z) - {{S}^{0}}$, где ${{S}^{0}}$ – удельная электропроводность вмещающей среды. Она может быть определена косвенными методами, например, через среднее кажущееся сопротивление.

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

где индексы означают: $p$ – точку наблюдения; $s$ –точку источника; $\nu $ – точку в неоднородности; $\Delta S$ – квадратная матрица электропроводности, имеющая вид: где: подматрицы $\Delta {{\sigma }_{j}} = \left( {\begin{array}{*{20}{c}} {\Delta \sigma _{{11}}^{j}}&{\Delta \sigma _{{12}}^{j}}&{\Delta \sigma _{{13}}^{j}} \\ {\Delta \sigma _{{21}}^{j}}&{\Delta \sigma _{{22}}^{j}}&{\Delta \sigma _{{23}}^{j}} \\ {\Delta \sigma _{{31}}^{j}}&{\Delta \sigma _{{32}}^{j}}&{\Delta \sigma _{{33}}^{j}} \end{array}} \right)$ есть тензоры избыточной электропроводности $j$-го элемента объема всей неоднородности $V$, состоящей из $N$элементов; ${\mathbf{X}}_{p}^{s}$ – первичное поле от источника в точке наблюдения; ${{{\mathbf{X}}}_{v}}$ – поле в неоднородности; ${{{\mathbf{X}}}^{p}}$ – поле в точке наблюдения (как разность потенциалов); $G_{v}^{p}$ – передаточная матрица от неоднородности в точку наблюдения; размерности $K \times 3N$, $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$ × $ \times \,\,{{([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}}}$.

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

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

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

получим:

Размерности матриц, входящих в последнюю систему уравнений, следующие: ${{[Z]}_{{m \times n}}}$_, ${{[X]}_{{m \times n}}}$, где $m$ – количество приемников с учетом размерности измеренного вектора в точке наблюдения; $n$ – количество источников с учетом размерности вектора стороннего тока в точке источника; ${{[Y]}_{{k \times n}}}$, ${{[\Delta S]}_{{k \times k}}}$, $[G_{v}^{p}] = m \times k$, ${{[G_{v}^{p}]}^{T}}_{{k \times m}}$, где $k$ – количество неоднородностей с учетом их размерности. Тогда

Перепишем (7) в виде $\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}}}$ имеют размеры 3 × 1, $M$ – количество источников. Для разрешимости каждого из этих уравнений необходимо минимум 3 источника.

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

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

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

Для анализа данной системы наблюдения рассмотрим задачу об определении анизотропии электропроводности в модели “колодец”. Данный объект представляет вертикально внедренное в нижнее однородное изотропное полупространство с удельной электропроводностью 1/30 См/м, анизотропное проводящее тело с тензором электропроводности:

Размеры этого тела 1 м × 1 м × 7 м, т.е. тело простирается на 7 м вниз.

Рассмотрим пример с симметричной системой наблюдения. Симметричной в том смысле, что источники, “освещающие” неоднородность, расположены симметрично относительно центра “колодца”. Система наблюдения и объект представлены на рис. 1.

Результаты моделирования системы наблюдения по указанному выше алгоритму приведены ниже.

Среднеквадратическое отклонение рассчитанного поля по данным решения обратной задачи от исходного поля для каждого источника, состоящего из пары заземленных электродов (в %): 7.437e–08, 1.569e–07, 1.264e–06, 1.512e–07. Величина среднеквадратического отклонения по электропроводности (в %): 51.126. Таким образом, при такой системе наблюдения точность решения обратной задачи относительно тензора анизотропии электропроводности оказывается невысокой – около 51%.

Рассмотрим систему наблюдения, асимметричную относительно изучаемого объекта, которая представлена на рис. 2.

Среднеквадратическое отклонение рассчитанного поля по данным решения обратной задачи от исходного поля для каждого источника, состоящего из пары заземленных электродов (в %): 8.625e–07, 8.338e–07, 2.741e–05, 6.315e–07. Величина среднеквадратического отклонения по электропроводности (в %): 2.152. Как следует из результатов моделирования, тензор удельной электропроводности восстанавливается с удовлетворительной точностью – примерно 2%.

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

ВЫВОДЫ

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

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

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

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