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

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

М. И. Шимелевич 1*, Е. А. Родионов 1**, И. Е. Оборнев 2***, Е. А. Оборнев 1****

1 Российский государственный геологоразведочный университет им. С. Орджоникидзе
г. Москва, Россия

2 Научно-исследовательский институт ядерной физики имени Д.В. Скобельцына МГУ им. М.В. Ломоносова
г. Москва, Россия

* E-mail: shimelevich-m@yandex.ru
** E-mail: evgeny_980@list.ru
*** E-mail: o_ivano@mail.ru
**** E-mail: obornevea@mail.ru

Поступила в редакцию 26.02.2022
После доработки 15.04.2022
Принята к публикации 15.04.2022

Полный текст (PDF)

Аннотация

В работе представлен аппроксимационный нейросетевой (АНС) метод решения обратной задачи геоэлектрики в кусочно-постоянных классах сред с последующим уточнением методом случайного поиска. Рассматривается условно-корректная обратная задача магнитотеллурического зондирования (МТЗ) на регуляризованной сетке параметризации. Число искомых параметров такой задачи может составлять $\sim {\kern 1pt} n \times {\text{1}}{{0}^{{\text{3}}}}$. Работа алгоритма иллюстрируется на примере инверсии натурных площадных наблюдений. Для полученного решения рассчитывались апостериорные характеристики его неоднозначности. Предлагаемые в работе методы реализованы с привлечением ресурсов суперкомпьютерных кластеров и технологии параллельных вычислений.

Ключевые слова: МТЗ, обратная задача, геофизика, нейронные сети, метод случайного поиска, натурные данные.

ВВЕДЕНИЕ

В работе представлен АНС метод решения обратной задачи геоэлектрики в кусочно-постоянных классах сред с последующим уточнением методом случайного поиска. АНС метод основан на аппроксимации обратного оператора при помощи нейросевых аппроксимационных конструкций. Рассматривается условно-корректная обратная задача магнитотеллурического зондирования (МТЗ) на регуляризованной сетке параметризации, блоки которой увеличиваются с глубиной. Критерием подбора числа ячеек сетки параметризации является “обучаемость” нейронной сети – на мелкой сетке параметризации сеть обучается плохо. Число параметров такой задачи может составлять $\sim {\kern 1pt} n \times {\text{1}}{{0}^{{\text{3}}}}$. Работа алгоритма иллюстрируется на примере инверсии натурных площадных наблюдений. Для полученного решения рассчитываются апостериорные характеристики его неоднозначности (погрешности). Предлагаемые в работе методы реализованы с привлечением ресурсов суперкомпьютерных кластеров и технологии параллельных вычислений.

ПОСТАНОВКА ЗАДАЧИ

Обратная задача МТЗ в заданном конечно-параметрическом классе сред сводится к решению операторного уравнения I рода вида [Шимелевич, Оборнев, 2009; Шимелевич и др., 2017]

(1)
${{A}_{N}}{\mathbf{s}} = {\mathbf{e}},\,\,\,\,{\mathbf{s}} \in {{S}_{N}} \subset {{R}^{N}},\,\,\,\,{\mathbf{e}} \in {{R}^{M}},\,\,\,\,M \geqslant N;$
где: ${\mathbf{s}} = ({{s}^{1}},...,{{s}^{N}})$ – искомый вектор параметров среды, каждый из которых может изменяться в заданном диапазоне значений $[s_{{min}}^{{}},s_{{min}}^{{}} + {{D}_{s}}]$; ${\mathbf{e}} = ({{e}^{1}},...,{{e}^{M}})$ – вектор входных данных; ${{A}_{N}}$ – непрерывный оператор численного решения прямой задачи в заданном конечно-параметрическом классе сред; ${{S}_{N}}$ – компактное множество векторов параметров, определяющих удельное сопротивление среды в исследуемой области $\Omega $ с учетом естественных априорных ограничений

${{S}_{N}}:\,\,\,[s_{{min}}^{{}} \leqslant {{s}_{n}} \leqslant s_{{min}}^{{}} + {{D}_{s}}],\,\,\,\,n = 1,...,N.$

РЕШЕНИЕ ЗАДАЧИ

Решение уравнения (1) есть некоторая неизвестная векторная функция от вектора данных ${\mathbf{e}}$:

${\mathbf{s}} = \psi ({\text{e}}),$

Предположим, что имеется множество решений (“банк решений”) ${{E}_{{bs}}}{\kern 1pt} :\,\, = \{ {{{\mathbf{e}}}^{i}},\;\;i = 1,...,{{Q}_{{bs}}}\} $ прямых задач, рассчитанных с помощью прямого оператора ${{A}_{N}}$ на некотором заданном множестве ${{S}_{{bs}}} \subset {{S}_{N}}$ векторов сред: ${{{\mathbf{e}}}^{i}} = {{A}_{N}}{{{\mathbf{s}}}^{i}}$, ${{{\mathbf{s}}}^{i}} \in {{S}_{{bs}}} \subset {{S}_{N}}$. Тогда для точных правых частей уравнения (1), принадлежащих ${{E}_{{bs}}}$, точные решения, т.е. точные значения функции $\psi $ в точках ${{{\mathbf{s}}}^{i}} \in {{S}_{{bs}}}$, известны:

${{{\mathbf{s}}}^{i}} = \psi ({{{\mathbf{e}}}^{i}}),\,\,\,\,i = 1,...,{{Q}_{{bs}}}.$

Аппроксимационный подход к решению уравнения (1) заключается в построении некоторой векторной функции заданного вида ${\mathbf{\Psi }}({\mathbf{\hat {a}}},{\mathbf{e}})$, определенной на множестве векторов данных ${\mathbf{e}} \in {{R}^{M}}$ и зависящей от матрицы коэффициентов ${\mathbf{\hat {a}}}$. В качестве функции ${\mathbf{\Psi }}({\mathbf{\hat {a}}},{\mathbf{e}})$ может быть задан, например, некоторый полином. Функцию ${\mathbf{\Psi }}({\mathbf{\hat {a}}},{\mathbf{e}})$ называют аппроксиматором обратного оператора для уравнения (1). Нахождение коэффициентов матрицы ${\mathbf{\hat {a}}}$ сводится к решению оптимизационной задачи приближения функции $\psi $ функцией $\Psi $ на множестве ${{E}_{{bs}}}$:

$\sum\limits_{i = 1}^{{{Q}_{{bs}}}} {\left\| {\Psi ({\mathbf{\hat {a},}}{{{\mathbf{e}}}^{i}}) - \psi ({{{\mathbf{e}}}^{i}})} \right\|_{{}}^{2}} \to \mathop {min}\limits_{{{\hat {a}}}} ,\,\,\,\,{{{\mathbf{e}}}^{i}} \in {{E}_{{bs}}}.$

Основным преимуществом данного метода является тот факт, что приближенное решение уравнения (1) может быть легко и быстро вычислено с помощью аппроксиматора $\Psi $ для любого вектора входных данных ${\mathbf{e}} \in {{R}^{M}}$. Построенный аппроксиматор тестируется на выборке известных решений, не принадлежащей множеству ${{E}_{{bs}}}$ (экзаменационной выборке), и затем вычисляется его собственная (интерполяционная) ошибка.

При АНС подходе к решению уравнения (1) в качестве функции $\Psi $ применяются нейросетевые полиномы или нейронные сети. Простейшей нейронной сетью является трехслойный персептрон, который описывается следующей формулой:

$\varphi ({\mathbf{V}},{\mathbf{W}},{\mathbf{x}}) = \sum\limits_{l = 1}^L {{{v}_{l}}} {\kern 1pt} \lambda \left( {\sum\limits_{m = 1}^M {{{w}_{{l,m}}}x_{m}^{{}}} } \right),$
где $\lambda $ – функция активации – заданная нелинейная, ограниченная, монотонно возрастающая, непрерывная функция, например, $\lambda (x) = {1 \mathord{\left/ {\vphantom {1 {(1 + \exp ( - x))}}} \right. \kern-0em} {(1 + \exp ( - x))}}$ (сигмоид). Вектор ${\mathbf{V}} = ({{v}_{1}},...,{{v}_{L}})$ и матрица ${\mathbf{W}} = \left\{ {{{w}_{{l,m}}}} \right\}_{{l,m = 1}}^{{L,M}}$ определяют матрицу ${\mathbf{\hat {a}}}$ нейросетевого аппроксиматора $\Psi ({\mathbf{\hat {a}}},{\mathbf{x}})$, а параметр L (количество нейронов на скрытом слое нейронной сети) определяет его сложность. На практике используются также нейросети, состоящие из пяти и более слоев. Использование нейронных сетей в задачах общей теории аппроксимации функций опирается на теоремы Колмогорова, Вейерштрасса–Стоуна, Цыбенко.

Многослойные нейросети – многослойные персептроны ($MLP$-сети), являются полносвязными сетями, в которых каждый элемент последующего слоя связан со всеми элементами предыдущего слоя, причем каждая связь имеет свой персональный весовой коэффициент. Последнее время в задачах анализа изображений широко применяются Глубокие нейронные сети (ГНС или DNN – Deep Neural Networks), состоящие из большого количества скрытых иерархических слоев с различными функциями нелинейных преобразований данных. Среди ГНС особое место занимают сверточные сети [Yoshua, 2009]. В сверточной нейронной сети каждый элемент последующего слоя данных связан с группой элементов (окном) предыдущего слоя (начиная с первого – входного слоя) посредством фильтра – операцией свертки; окно, образуя ядро свертки, скользит по слою, в самом начале – непосредственно по первому входному слою, то есть анализируемому цифровому массиву входных данных. После многократных прохождений указанных слоев преобразований данные передаются на обычную полносвязную нейронную сеть. В задачах анализа цветных изображений первым входным слоем сети является трехмерный массив данных.

В работах [Шимелевич, Оборнев, 2009; Шимелевич и др., 2017; Shimelevich et al., 2018] показана эффективность использования локальных НС аппроксиматоров при решении обратных задач геофизики, которые устанавливают количественную связь некоторой группой данных (окном данных входного слоя) с определенной группой параметров (окном параметров) среды. При этом производится сжатие данных (частичная инверсия на основе неполных данных). При этом, если измеряется несколько характеристик поля на сетке частот, то входной слой является “тензорным”, так как включает набор трехмерных массивов на каждой частоте для каждой характеристки. Глубина тензорного входного слоя кратна числу частот и числу используемых характеристик поля. После указанных преобразований данные передаются на полносвязную $MLP$-сеть.

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

ПОСТРОЕНИЕ НЕЙОСЕТЕВОГО АППРОКСИМАТОРА

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

Рис. 1.

Вертикальный разрез модели M-3D. Тонкими линиями обозначена сетка прямой задача, жирными – сетка параметризации.

Параметрами модели M-3D являются десятичные логарифмы удельных сопротивлений $\lg \rho _{n}^{{}}$ в центрах блоков параметризации. Как видно из рис. 1, детальность параметризации среды с глубиной уменьшается в силу затухания электромагнитного поля, так как прямая задача МТЗ описывается системой уравнений Максвелла в квазистационарном приближении и отклик среды от более глубоких фрагментов среды слабее, чем от лежащих ближе к поверхностности Земли. Критерием при определении размеров блоков в данной работе была выбрана обучаемость нейронной сети с приемлемой ошибкой на экзаменационной выборке (до 12%): для каждого яруса сетки параметризации подбиралось максимальное число блоков, при котором нейронная сеть обучается с приемлемой ошибкой. Общее количество параметров модели составило 1225. Данная сетка параметризации, покрывающая исследуемую область $\Omega $, размеры области диапазон $[s_{{min}}^{{}},s_{{min}}^{{}} + {{D}_{s}}]$ изменения искомых параметров среды задают класс сред $G_{1}^{{3D}}$ с кусочно-постоянным распределением уд. сопротивления, в котором решается обратная задача (1).

Для модельного класса сред $G_{1}^{{3D}}$ рассчитывался банк решений прямых задач, при этом логарифм удельного сопротивления параметризованных блоков имел случайное равномерное распределение в диапазоне от 0 до 4. В качестве численного прямого оператора использовался модуль из известной программы R. Mackie mtd3fwd, находящейся в свободном доступе.

В рамках банка решений было рассчитано 50 000 прямых задач на сетке из 14 периодов в диапазоне от 0.001 с до 21.5 с. Время расчета одной прямой задачи на данной сетке периодов на одном CPU персонального компьютера составляет примерно 30 минут. Вследствие высокой вычислительной сложности банк решений рассчитывался с привлечением ресурсов суперкомпьютерных кластеров и технологий параллельных вычислений. При расчете прямой задачи на каждом из периодов рассчитывался вектор данных, состоящий из действительных и мнимых частей четырех компонент тензора импеданса в узлах пространственной XY-сетки, состоящей из 900 узлов.

Для данной модели были построены локальные пятислойные нейросетевые аппроксиматоры с количеством нейронов на трех скрытых слоях 32, 16 и 8. Аппроксиматоры тестировались на независимой экзаменационной выборке известных решений обратных задач, отсутствующих в банке решений. Для l-го слоя сетки вычислялась собственная ошибка аппроксиматоров (табл. 1) как средняя ошибка по параметрам слоя и по выборке, деленная на диапазон изменения параметров:

$\varepsilon _{\Psi }^{l} = \frac{1}{{4{{S}_{T}}{{J}_{l}}}}\sum\limits_{i = 1}^{{{S}_{E}}} {\sum\limits_{j = 1}^{{{J}_{l}}} {\left| {\Psi ({\mathbf{\hat {a}}},{{{\mathbf{e}}}^{i}}) - \lg \rho _{j}^{i}} \right|} } ,$
где: $l = 1,...,5,$ – номер слоя; ${{J}_{l}}$ – количество параметров на l-м слое; ${{S}_{E}}$ – объем экзаменационной выборки.

Таблица 1.  

Собственные ошибки нейросетевых аппроксиматоров для класса сред M-3D

№ слоя 1 2 3 4 5
Средняя ошибка по экзаменационной выборке, % 1.16 2.89 4.3 7.17 10.8

Как видно из табл. 1, ошибка аппроксиматоров увеличивается с глубиной, но в целом не превышает 11%. Обученный аппроксиматор позволяет решать обратные 3D-задачи, сводящиеся к уравнению (1) в классе сред $G_{1}^{{3D}}$ за время в пределах первых десятков секунд и, таким образом, проводить интерпретацию измеренных данных в режиме реального времени.

Построенные аппроксиматоры были успешно оттестированы также на модельных примерах сред из класса $G_{1}^{{3D}}$ [Shimelevich et al., 2018].

МЕТОД УТОЧНЕНИЯ РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ

Следует отметить, что АНС метод является универсальным для получения приближенного решения обратной задачи в заданном классе сред, но при этом он не минимизирует невязку. Получаемое приближенное нейросетевое решение, при необходимости, можно уточнять различными методами. Например, в работе [Шимелевич и др., 2017] предложен аппроксимационно-итерационный нейросетевой метод. Однако в случае высокой размерности сетки прямой задачи данный метод является очень времяемким с вычислительной точки зрения. Менее затратным является метод случайного поиска, который заключается в поиске группы решений в окрестности найденного, с невязкой синтеза меньшей, чем исходная. При этом по полученному распределению логарифмов удельных сопротивлений стандартным методом кластеризации k-means определяются геометрические границы схожести параметров. Случайный поиск осуществляется с учетом этих границ: в окрестности найденного решения насчитывается множество новых решений, причем изменение параметров относительно исходного вектора происходит на уровне групп параметров – кластеров. Оптимальное число кластеров находилось при помощи метода псевдо-F-статистики [Любушин, 2011]. Предположим, что в окрестности первого приближения найдено K векторов решений ${{{\mathbf{s}}}^{l}} = (s_{1}^{l},s_{2}^{l},...,s_{n}^{l}),$ $l = 1,...,K,$ с невязками ${{\delta }_{l}}$ меньшими, чем исходная, и отсортированными по возрастанию. Тогда новое множество векторов решений вычисляется по формулам:

$\begin{gathered} (s_{1}^{{'l}},s_{2}^{{'l}},...,s_{n}^{{'l}}) = \frac{1}{K}\sum\limits_{k = 1}^K {\delta (k)\sum\limits_{l = 1}^k {\frac{1}{{{{\delta }_{l}}}}} } (s_{1}^{l},s_{2}^{l},...,s_{n}^{l}), \\ \delta (k) = \sum\limits_{l = 1}^k {{{\delta }_{l}},} \,\,\,\,l = 1,...,K. \\ \end{gathered} $
Для каждого из них снова вычисляется невязка синтеза, а за приближенное решение принимается вектор, доставляющий наименьшую невязку.

В вычислительной практике наряду с классическим расстоянием в нормированном пространстве X, определяемым нормой разности $\left\| {x - x{\kern 1pt} '} \right\|$, в качестве меры близости элементов $x,\,x{\kern 1pt} ' \in X$ бывает целесообразно оценивать относительное расстояние $d_{X}^{r}(x,x{\kern 1pt} ')$ между элементами. Например, часто используется функционал относительного расстояния, между $x$ и $x{\kern 1pt} '$ вида [Беклемишев, 1983; Годунов и др., 1988]:

(2)
$d_{X}^{r}(x{\kern 1pt} ',x) = \frac{{\left\| {x{\kern 1pt} '\,\, - x} \right\|}}{{\left\| x \right\|}},\,\,\,\,x \ne 0.$

Определяемое функционалом $d_{X}^{r}(x{\kern 1pt} ',x)$ относительное расстояние типа (2) не является метрикой в нормированном пространстве X, однако по физическому смыслу эта величина является количественной мерой близости элементов $x$ и $x{\kern 1pt} '$.

Обозначим через ${{\tilde {e}}_{i}}(j,t)\,$ и $e_{i}^{'}(j,t)$ комплексные значения измеренной и синтезированной j-й электромагнитной характеристики (одной из четырех компонент тензора импеданса: Zxx, Zxy, Zyx, Zyy), соответственно, в i-ой точке на t-ом периоде, $i = 1,...,n,$$j = 1,...,4,$ $t = 1,...,Nt.$ Здесь n – количество точек пространственной XY-сетки, Nt – число периодов. Квадраты дифференцированных относительных площадных невязок вычисляются по формулам:

(3)
$\begin{gathered} \delta _{{j,t}}^{2} = \frac{{\sum\limits_{i = 1}^n {{{{\left| {{{{\tilde {e}}}_{i}}(j,t) - e_{i}^{'}(j,t)} \right|}}^{2}}} }}{{\sum\limits_{i = 1}^n {{{{\left| {{{{\tilde {e}}}_{i}}(j,t)} \right|}}^{2}}} }}, \\ j = 1,...,4,\,\,\,\,t = 1,...,Nt. \\ \end{gathered} $

Введем теперь весовые коэффициенты для четырех электромагнитных характеристик:

(4)
${{w}_{j}} = \sum\limits_{t = 1}^{Nt} {\sqrt {\frac{1}{n}\sum\limits_{i = 1}^n {{{{\left| {\tilde {e}(j,t)} \right|}}^{2}}} } } ,\,\,\,\,j = 1,...,4.$

Тогда обобщенная площадная невязка с весовыми коэффициентами вычисляется по формуле:

(5)
$\delta = {{\left( {\sum\limits_{j = 1}^4 {{{w}_{j}}} } \right)}^{{ - 1}}}\sum\limits_{j = 1}^4 {{{w}_{j}}} {{\left( {\frac{1}{{Nt}}\sum\limits_{t = 1}^{Nt} {\delta _{{j,t}}^{2}} } \right)}^{{1/2}}}.$

Весовые коэффициенты (4) определяют значимость каждой из компонент тензора импеданса в соответствии с уровнем измеренных характеристик по всему планшету измерений. Величина, определяемая формулой (5), не является нормой, а представляет собой функционал относительного расстояния $d({\mathbf{\tilde {e}}},{\mathbf{e}}{\kern 1pt} ')$ с весовыми множителями в пространстве F, где ${\mathbf{\tilde {e}}},{\mathbf{e}}{\kern 1pt} '$ – измеренный и синтезированный векторы данных, соответственно.

АПОСТЕРИОРНЫЕ ОЦЕНКИ НЕОДНОЗНАЧНОСТИ РЕШЕНИЯ

Для условно-корректной обратной задачи может быть получена апостериорная оценка степени неоднозначности (погрешности) получаемого приближенного решения [Шимелевич и др., 2018]. Положим ${\mathbf{\tilde {e}}} = {\mathbf{e}} + \Delta {\mathbf{e}}$, где $\Delta {\mathbf{e}}$ – вектор погрешности правой части в уравнении (1) такой, что ${{\left\| {\Delta {\mathbf{e}}} \right\|}_{F}} = {{\delta }_{0}}$. Обозначим через ${{{\mathbf{s}}}^{0}}$ приближенное решение, удовлетворяющее уравнению

(6)
${{{A}_{N}}{\text{s}} = {{\tilde {e}}},\,\,\,\,{\text{s}} \in {{G}_{1}}} \subset G,\,\,\,\,{{\tilde {e}}} \in F,$
с невязкой, не превышающей величину ${{\delta }_{0}}$ погрешности в правой части уравнения (6):

$\left\| {{{A}_{N}}{{{\mathbf{s}}}^{0}} - {\mathbf{\tilde {e}}}} \right\| \leqslant {{\delta }_{0}}.$

Максимально возможное отклонение любого другого приближенного решения ${{{\mathbf{s}}}^{1}}$ такого, что $\left\| {{{A}_{N}}{{{\mathbf{s}}}^{1}} - {\mathbf{\tilde {e}}}} \right\| \leqslant {{\delta }_{0}}$ (т.е. ${{\delta }_{0}}$ – эквивалентного решения), от решения ${{{\mathbf{s}}}^{0}}$ не может превышать величины $\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},{{\delta }_{0}})$:

$\left\| {{{s}^{1}} - {{s}^{0}}} \right\| \leqslant \beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},{{\delta }_{0}}),$
где величина $\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},{{\delta }_{0}})$ определяется выражением:

(7)
$\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},{{\delta }_{0}}) = \mathop {\max }\limits_{s \in {{S}_{N}}} \left\{ {\left\| {{\mathbf{s}} - {{{\mathbf{s}}}^{0}}} \right\|\,\,:\,\,\left\| {{{A}_{N}}{\mathbf{s}} - {\mathbf{\tilde {e}}}} \right\| \leqslant {{\delta }_{0}}} \right\}.$

Если уровень ${{\delta }_{0}}$ погрешности правой части ${\mathbf{\tilde {e}}}$ известен, а ${{{\mathbf{s}}}^{{Ex}}}$ – точное решение уравнения (1), то для ${{{\mathbf{s}}}^{{Ex}}}$ выполняется условие $\left\| {{{A}_{N}}{{{\mathbf{s}}}^{{Ex}}} - {\mathbf{\tilde {e}}}} \right\| \leqslant {{\delta }_{0}}$, и может быть получена апостериорная оценка погрешности $\varepsilon ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},{{\delta }_{0}})$ найденного решения ${{{\mathbf{s}}}^{0}}$ уравнения (6) для фиксированной правой части ${\mathbf{\tilde {e}}}$ [Гончарский, Черепащук, Ягола, 1978]:

(8)
$\varepsilon ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},{{\delta }_{0}}) = \left\| {{{{\mathbf{s}}}^{{Ex}}} - {{{\mathbf{s}}}^{0}}} \right\| \leqslant \beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},{{\delta }_{0}}).$

На практике в задачах геофизики уровень погрешности данных ${{\delta }_{0}}$, как правило, неизвестен, а известна фактическая невязка $\delta $ найденного приближенного решения ${{{\mathbf{s}}}^{0}} \in {{S}_{N}}$. В этом случае неравенство (8) при ${{\delta }_{0}} = \delta $ определяет степень неоднозначности приближенного решения ${{{\mathbf{s}}}^{0}}$, полученного для фиксированной правой части ${\mathbf{\tilde {e}}}$ с невязкой $\delta $:

$\begin{gathered} \left\| {{\mathbf{s}} - {{{\mathbf{s}}}^{0}}} \right\| \leqslant \beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},{{\delta }_{0}}),\,\,\,\,{\text{для всех}}~ \\ {\mathbf{s}} \in {{S}_{N}}\,\,:\,\,\,\left\| {{{A}_{N}}{\mathbf{s}} - {\mathbf{\tilde {e}}}} \right\| \leqslant \delta , \\ \end{gathered} $
а величина $\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},{{\delta }_{0}})$ называется апостериорной оценкой степени неоднозначности решения ${{{\mathbf{s}}}^{0}}$.

На практике возникает потребность оценить степень неоднозначности решения в некоторой подобласти ${{\Omega }_{1}} \subset \Omega $, определяющей свойства некоторой выделенной геологической структуры. Обозначим ${\mathbf{s}} = {{{\mathbf{s}}}^{0}} + \Delta {\mathbf{s}}$. Если компоненты вектора $\Delta {\mathbf{s}}$ отличны от нуля лишь на множестве индексов, отвечающем подобласти ${{\Omega }_{1}}$, то формула (7) определяет локальную апостериорную оценку неоднозначности полученного решения для этой подобласти.

Локальные апостериорные оценки неоднозначности решения можно найти численно методом Монте-Карло.

Перепишем (7) в эквивалентном виде, заменив ${{\delta }_{0}}$ на значение величины невязки синтеза $\delta $ полученного решения ${{s}^{0}}$:

(9)
$\begin{gathered} \beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},\delta ) = \mathop {\max }\limits_{} \left\{ {\left\| {\Delta {\mathbf{s}}} \right\|\,\,:\,\,\left\| {\Delta {\mathbf{e}}({{{\mathbf{s}}}^{0}},\Delta {\mathbf{s}})} \right\| \leqslant \delta ,} \right. \\ \left. {\Delta {\mathbf{e}}({{{\mathbf{s}}}^{0}},\Delta {\mathbf{s}}) = {{A}_{N}}({{s}^{0}} + \,\,\Delta {\mathbf{s}}) - {\mathbf{\tilde {e}}},\,\,{{{\mathbf{s}}}^{0}} + \,\,\Delta {\mathbf{s}} \in {{S}_{N}}} \right\}. \\ \end{gathered} $

Из (9) следует, что величина $\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},\delta )$ определяется парой скалярных величин $\left\| {\Delta {\mathbf{s}}} \right\|$ и $\left\| {\Delta {\mathbf{e}}({{{\mathbf{s}}}^{0}},\Delta {\mathbf{s}})} \right\|$, принимающих различные значения при всех возможных значениях вектора $\Delta {\mathbf{s}}$. Всевозможные пары $\left\| {\Delta {\mathbf{s}}} \right\|$ и $\left\| {\Delta {\mathbf{e}}({{{\mathbf{s}}}^{0}},\Delta {\mathbf{s}})} \right\|$ определяют множество точек на декартовой плоскости:

$\begin{gathered} \Lambda = \left\{ {\left( {~\left\| {\Delta {\mathbf{s}}} \right\|,\,\left\| {\Delta {\mathbf{e}}({{{\mathbf{s}}}^{0}},\Delta {\mathbf{s}})} \right\|} \right)\,\,:\,\,\Delta {\mathbf{e}}({{{\mathbf{s}}}^{0}},\Delta {\mathbf{s}})} \right. = \\ \left. { = {{A}_{N}}({{{\mathbf{s}}}^{0}} + \Delta {\mathbf{s}}) - {\mathbf{\tilde {e}}},\,\,{{{\mathbf{s}}}^{0}} + \,\Delta {\mathbf{s}} \in {{G}_{1}}[\Omega ]} \right\}. \\ \end{gathered} $

Если множество $\Lambda $ построено, то проблема отыскания величины $\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},\delta )$ сводится к нахождению точки с максимальной абсциссой на подмножестве $\Lambda $, удовлетворяющем условию $\left\| {\Delta {\mathbf{e}}({{{\mathbf{s}}}^{0}},\Delta {\mathbf{s}})} \right\| \leqslant \delta $. Абсцисса этой точки и будет величиной $\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},\delta )$.

Для численного нахождения $\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},\delta )$ с по-мощью прямого оператора во всем допустимом диапазоне насчитывается множество $\tilde {\Lambda } \subset \Lambda $, на котором приближенно находится апостериорная оценка неоднозначности решения. В работе [Шимелевич, 2020] доказана сходимость по вероятности этого алгоритма.

Характеристики степени неоднозначности для различных функционалов расстояния рассматривались, например, в работах [Канторович, 1962; Беклемишев, 1983; Новик, 1987; Годунов и др., 1988; Шимелевич, 1989] и ряде других работ.

РЕЗУЛЬТАТЫ

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

Рис. 2.

Схема расположения профилей магнитотеллурической съемки участка Северо-Пясинской площади.

Характеристиками МТЗ данной съемки являются действительные и мнимые части четырех компонент тензора импеданса в диапазоне периодов от 0.0046 до 631 с. Съемка представляет собой 10 профилей на расстоянии примерно 1 км друг от друга, расстояние между пикетами также примерно 1 км. Размер исследуемого участка $9.5 \times 10$ км. Эти данные прошли первичную обработку, заключающуюся в фильтрации и нормализации [Фельдман и др., 2008].

После получения первого нейросетевого приближения было проделано еще 5 итераций методом случайного поиска, в ходе которого для вычисления невязок применялась формула (5). На рис. 3 показано, как уменьшается невязка при увеличении номера итерации. На каждой итерации объем выборки случайного поиска составил около 2000 решений прямых задач.

Рис. 3.

График уменьшения невязки $\delta $ по итерациям.

Рис. 4.

Иллюстрация численного нахождения апостериорной оценки неоднозначности решения на шестой итерации для пятого слоя.

Рис. 5.

Горизонтальные сечения решения на шестой итерации (а) и решения, отвечающего точке множества $\tilde {\Lambda }$ с минимальной ординатой ($24.72\% $), для первого слоя.

На последней итерации невязка $\delta $ составила 24.70%. На рис. 6а–рис. 10а представлены горизонтальные сечения пяти слоев решения, полученного на шестой итерации. Цветом обозначен десятичный логарифм удельного сопротивления.

Рис. 6.

Горизонтальные сечения решения на шестой итерации (а) и решения, доставляющего оценку $\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},\delta )$, для второго слоя.

Рис. 7.

Горизонтальные сечения решения на шестой итерации (а) и решения, доставляющего оценку $\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},\delta )$, для третьего слоя.

Рис. 8.

Горизонтальные сечения решения на шестой итерации (а) и решения, доставляющего оценку $\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},\delta )$, для четвертого слоя.

Рис. 9.

Горизонтальные сечения решения на шестой итерации (а) и решения, доставляющего оценку $\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},\delta )$, для пятого слоя.

Для решения на шестой итерации, соответствующего наблюдениям Северо-Пясинской площади, были рассчитаны послойные апостериорные оценки его неоднозначности для пяти слоев, из которых состоит модель M-3D. При этом в пространстве параметров была выбрана следующая норма:

$\left\| {\Delta {\mathbf{s}}} \right\| = \frac{1}{4}\max \{ {{\left( {\Delta {\mathbf{s}}} \right)}_{n}},\,\,\,\,n = 1,...,N\} ,$
где множитель $\frac{1}{4}$ выбран в соответствии с принятым диапазоном изменения логарифма удельного сопротивления среды. Вместо величины $\left\| {\Delta {\mathbf{e}}({{{\mathbf{s}}}^{0}},\Delta {\mathbf{s}})} \right\| = \left\| {{\mathbf{\tilde {e}}} - {{A}_{N}}({{{\mathbf{s}}}^{0}} + \Delta {\mathbf{s}})} \right\|$ вычислялся функционал относительного расстояния $d({\mathbf{\tilde {e}}},{{A}_{N}}{\mathbf{s}})$, определяемый формулами (3)(5), в которых в качестве $e_{i}^{'}(j,t)$ выступают i-е компоненты комплексных значений вектора ${{A}_{N}}({{{\mathbf{s}}}^{0}} + \Delta {\mathbf{s}})$ j-й электромагнитной характеристики на t-м периоде. В этом случае апостериорная оценка неоднозначности полученного решения ${{{\mathbf{s}}}^{0}}$ будет определяться выражением:

$\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},\delta ) = \mathop {\max }\limits_{s \in {{G}_{1}}[\Omega ]} \left\{ {\left\| {{\mathbf{s}} - {{{\mathbf{s}}}^{0}}} \right\|\,\,:\,\,d({\mathbf{\tilde {e}}},{{A}_{N}}{\mathbf{s}}) \leqslant \delta } \right\}.$

Метод нахождения апостериорной оценки неоднозначности решения на шестой итерации для пятого слоя проиллюстрирован на рис. 4. Послойные апостериорные оценки неоднозначности приведены в табл. 2. Для первого слоя все построенное множество $\tilde {\Lambda }$, оказалось выше линии невязки $\delta = 0.247$, вследствие чего определить значение величины $\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},\delta )$ для первого слоя не удалось, и была определена “квазиоценка” для решения, отвечающего точке множества $\tilde {\Lambda }$ с минимальной ординатой ($24.72\% $). На рис. 5рис. 9 представлены горизонтальные сечения найденного на шестой итерации решения (а) и решения, доставляющие оценку $\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},\delta )$ для каждого из слоев модели M-3D (б), то есть наиболее удаленные $\delta $‑эквивалентные решение при полученном значении невязки синтеза $\delta $; цветом обозначен десятичный логарифм удельного сопротивления.

Таблица 2.  

Послойные апостериорные оценки степени неоднозначности решения, полученного на шестой итерации и соответствующие им невязки синтеза

№ слоя 1 2 3 4 5
Апостериорная оценка $\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},\delta = 0.247)$, % <0.5 12.39 3.15 4.15 14.36
Невязка $\delta $-эквивалентного решения, доставляющего оценку $\beta ({\mathbf{\tilde {e}}},{{{\mathbf{s}}}^{0}},\delta )$, % 24.69 24.69 24.67 24.51

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

Получены первые удовлетворительные результаты 3D-нейросетевой инверсии натурных данных с числом параметров порядка n × 103. Построенный нейросетевой аппроксиматор может быть многократно и единообразно применен для решения обратных 3D-задач МТЗ (аналогично использованию электронных 1D-палеток Ваньяна) с получением результатов инверсии в режиме реального времени. Первое нейросетевое приближение многократно уточнялось методом случайного поиска с применением алгоритма кластериации k-means. Получены апостериорные оценки неоднозначности для решения на последней итерации, определяющие возможный разброс $\delta $-эквивалентных решений при достигнутой невязке.

Список литературы

  1. Беклемишев Д.В. Дополнительные главы линейной алгебры. М.: Наука. 1983. 336 с.

  2. Годунов С.К., Антонов А.Г., Кирилюк О.П., Костин В.И. Гарантированная точность решения систем линейных уравнений в евклидовых пространствах. М: Наука. 1988.456 с.

  3. Гончарский А.В., Черепащук А.М., Ягола А.Г. Численные методы решения обратных задач астрофизики. М.: Наука. 1978. 36 с.

  4. Канторович Л.В. О некоторых новых подходах к вычислительным методам и обработке наблюдений // Сибирский математический журн. 1962. Т. Ш. № 5. С. 701–709.

  5. Любушин А.А. Кластерный анализ свойств низкочастотного микросейсмического шума // Физика Земли. 2011. № 6. С. 26–34.

  6. Новик О.Б. Математические вопросы сокращения числовой геофизической информации при поисках нефти и газа. Деп. в ВИЭМС 02.11.87. № 485-МГ.

  7. Шимелевич М.И. Методы решения обратных задач электромагнитных зондирований в двумерном приближении. Деп. в ВИЭМС 02.12.88 № 670-МГ-89.

  8. Шимелевич М.И. О методе расчета модуля непрерывности обратного оператора и его модификаций с приложением к нелинейным задачам геоэлектрики // Вычислительные методы и программирование. 2020. Т. 21. Вып. 4. С. 350–372.

  9. Шимелевич М.И., Оборнев Е.А. Аппроксимационный метод решения обратной задачи МТЗ с использованием нейронных сетей // Физика Земли. 2009. № 12. С. 22–38.

  10. Шимелевич М.И., Оборнев Е.А., Оборнев И.Е., Родионов Е.А. Аппроксимационный нейросетевой метод решения многомерных нелинейных обратных задач геофизики // Физика Земли. 2017. №4. С. 100–109.

  11. Шимелевич М.И., Оборнев Е.А., Оборнев И.Е., Родионов Е.А. Алгоритм решения обратной задачи геоэлектрики на основе нейросетевой аппроксимации // Сибирский журнал вычислительной математики. 2018. № 4. С. 437–452.

  12. Фельдман И.С., Окулесский Б.А., Сулейманов А.К., Николаева В.И., Кунчеров В.А., Чамо С.С. Электроразведка методом МТЗ в комплексе региональных нефтегазопоисковых работ в европейской части России // Записки Горного института. 2008. Т. 176. С. 125–131.

  13. Bengio Y. Learning Deep Architectures for AI // Foundations and Trends in Machine Learning. 2009. V. 2. № 1. P. 1–127.

Дополнительные материалы отсутствуют.