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

ВВЕДЕНИЕ

Интерпретацию экспериментальных данных нельзя считать завершенной без оценки качества полученных результатов; его принято характеризовать точностью приближенного решения обратной задачи, отвечающего некоторому критерию оптимальности. Для операторных уравнений проблема построения оценок точности – как априорных, так и апостериорных – рассматривалась, в том числе, в работах [Леонов, 2010; 2014]. Уязвимым местом теоретических оценок в плане практического применения является присутствующая в них неизвестная константа, для установления приближенного значения которой одних математических приемов может оказаться недостаточно и приходится прибегать к имитационному моделированию [Барашков, 1993]. В теории интерпретации гравитационных и магнитных аномалий критерии надежности и точности приближенных решений обратных задач являются неотъемлемой частью физико-математической интерпретационной модели [Страхов, 1995]. На практике же большинство известных алгоритмов решения обратных задач гравиразведки не предусматривают оценку качества своих результатов, что объясняется ограниченными возможностями общепринятых математических форм представления результатов интерпретации и чрезмерной абсолютизацией свойства оптимальности.

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

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

где $S{\text{*}}$ и $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} $ суть построенное приближенное и неизвестное точное решения обратной задачи. В практике решения обратных задачах гравиразведки оценка (1) используется преимущественно для характеристики работы того или иного алгоритма на модельных примерах. Принципиально отличными от (1) являются непосредственные оценки истинного решения обратной задачи, не привязанные ни к одному конкретному допустимому варианту интерпретации. Примером подобных оценок являются минимально и максимально возможные значения эффективной массы источников гравитационного поля, которые акад. Л.В. Канторович предлагал рассматривать как новую самостоятельную форму представления решения линейной обратной задачи гравиразведки [Канторович, 1962]. В обратной задаче гравиразведки рудного типа к альтернативному типу оценок можно отнести минимальную и максимальную пространственные области ${{D}_{1}}$ и ${{D}_{2}}$, обеспечивающие неулучшаемые включения для носителя $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} $ источников аномалии. Пара $({{D}_{1}},{{D}_{2}})$ была предложена в работе [Балк, 1980] как особая форма представления результатов интерпретации в обратной задаче для одиночного рудного тела. В настоящей статье приводится конструктивный и полностью формализованный метод построения раздельных двухсторонних оценок типа (2) для локальных (парциальных) носителей в описании многосвязного распределения аномалиеобразующих масс при ограничениях, охватывающих едва ли не все типы априорной информации, встречающихся в известных алгоритмах решения обратной задачи гравиразведки. Последнее оказалось возможным благодаря геофизически ориентированному (монтажному) методу поиска единичных приближенных решений обратной задачи в конечноэлементных классах источников поля, используемому при построении состоятельного подмножества допустимых вариантов интерпретации.

ОСНОВНЫЕ ОТЛИЧИТЕЛЬНЫЕ ОСОБЕННОСТИ ОЦЕНОК (1) И (2)

Отметим некоторые принципиальные с методологической точки зрения особенности оценок двух типов. Начнем с оценок типа (1). Практическая ценность их – что вполне естественно – во многом будет определяться значимостью и информативностью математического объекта, который обслуживают эти оценки, то есть приближенного решения $S{\text{*}}$ обратной задачи. Здесь необходимо отметить следующее. Множество (в нашем случае, множество $Q$ допустимых решений обратной задачи) невозможно сколь-нибудь полно охарактеризовать одним из его элементов, какими бы замечательными свойствами он ни обладал. Соответственно, информативность единичного решения обратной задачи может значительно уступать информативности всего множества $Q$ – результаты интерпретации в терминах отдельного решения обратной задачи не позволяют претендовать на полноту извлечения информации из геофизических данных. Статус лучшего, который обычно приписывается приближенному решению $S{\text{*}}$, ситуацию не спасает, не говоря уже о том, что не всегда декларируемая оптимальность действительно отвечает за точность решения (взять, хотя бы, “оптимальность по минимуму достигнутой невязки”). Как отмечал акад. Я.З. Цыпкин [Цыпкин, 1992], в отношении методов оптимальной обработки данных оправдывается афоризм Дж. Тьюки “оптимальность становится опасной, если ее принимать слишком всерьез”. Фактическое качество приближенного решения $S{\text{*}}$, по какому бы критерию оптимальности оно ни было выбрано, является случайным настолько, насколько случайна выборка помех в измерениях поля. Это не позволяет обеспечить монотонную зависимость качества приближенного решения от объема априорной информации. Как следствие из сказанного выше напрашивается вывод: не самый лучший прием использовать в конструкции оценок математический объект со столь нестабильными свойствами. К сожалению, базовое для теории некорректных обратных задач положение о (не реализуемой на практике) возможности алгоритма получать решения с наперед заданной точностью во многом скрывают отмеченные недостатки.

И это еще не все. Оценки (1) дают опосредованное представление об истинном решении $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} $. Если точность приближенного решения не высока и не позволяет отождествить его с истинным, возникает дополнительная задача восстановления свойств истинного решения по паре ($S{\text{*}}$, $\Delta $). Сделать это не всегда просто. Взять, к примеру, задачу приближенного восстановления геометрического места точек, принадлежащих носителю $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} $ в случае, когда обратная задача была поставлена в классе многогранников. К тому же, и это главное, свертывание информации о носителе $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} $ в пару $(S{\text{*}},\Delta )$ приводит к ее потере, поскольку неравенству (1) помимо допустимых удовлетворяют и недопустимые решения – если оценку (1) “развернуть” в оценку (2), то последняя может оказаться значительно грубее, построенной непосредственно по множеству $Q$ приближенных решений решений обратной задачи, отвечающих всем априорным данным и предпосылкам. Поясним это. Пусть скалярная величина $x$ априори принадлежит интервалу $Q = \left[ {0;1} \right]$ и, исходя из каких-то соображений, в качестве приближенного выбрано значение $x* = 0.2$. Тогда из неулучшаемой оценки $\left| {x - 0.2} \right| \leqslant 0.8$ следует, что $x \in \left[ { - 0.6;1} \right]$.

Обратимся теперь к оценкам (2). Прежде всего, они свободны от субъективных моментов, связанных с назначением критерия оптимальности решения $S{\text{*}}$. Свойство оптимальности, равно как и отдельные решения, априори наделенные преимуществами перед другими, вообще не фигурируют в оценках (2). В обратных задачах гравиразведки, где в силу дискретности измерений поля используются конечнопараметрические модельные классы $M$, множество $Q$ суть ограниченное множество из R$^{n}$. Его “размеры” суть степень распространения приближенной эквивалентности. Смещение акцента с единичного решения $S{\text{*}}$ на множество допустимости $Q$ позволяет поднять интерпретацию на уровень анализа объективных возможностей геофизического метода, а не отдельного взятого алгоритма решения обратной задачи.

Оценки (2) напоминают о неустоявшейся терминологии в отношении задач локализации и детализации геологического объекта по геофизическим данным. Какую конкретно из этих двух задач удалось решить, интерпретатор зачастую определяет по предполагаемому значению достигнутой точности оценивания параметров модели источников поля. Если (не претендуя на общность введенных понятий, а исходя из общепринятого смысла этих слов) под локализацией геологического объекта понимать построение (по возможности минимальной) области, в пределах которой сосредоточены все (или основные) аномальные массы, а под детализацией – построение (по возможности максимального) фрагмента носителя $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} $, то оценки ${{D}_{1}}$ и ${{D}_{2}}$, как раз, и будут обслуживать эти задачи.

Важно подчеркнуть, что каждая из оценок (1) и (2) не перекрывает возможности другой и несет свою информацию об изучаемом объекте, и потому их надо рассматривать не как конкурирующие, а скорее как дополняющие друг друга. Несомненное достоинство оценок (2) состоит в достоверности доставляемой ими информации о пространственном положении источников аномалии (разумеется, при условии адекватности всех априорных посылок), что не свойственно оценкам типа (1). Однако при довольно скудном объеме априорных данных такая информация может оказаться малозначимой, а то и вовсе нулевой. Последнее произойдет тогда, когда область ${{D}_{2}}$ совпадает с известным фрагментом истинного носителя $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} $, а область ${{D}_{1}}$ совпадает с областью, априори содержащей источники аномалии. В той же самой ситуации при использовании порожденной мерой Лебега $\mu $ метрики Штейнхауса [Marczewski, Steinhaus, 1958] :

полезную информацию можно извлечь уже из оценки $\rho (S*,S) \leqslant \Delta $. Так, можно оценить константу $\alpha \in $[0,1]: $\mu (S{\text{*}} \cap \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} ) = \alpha \mu (\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} )$, иначе говоря, оценить меру фрагмента, общего для носителей $S{\text{*}}$ и $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} $, по отношению к мере последнего, хотя при этом нельзя уточнить, о каком конкретно фрагменте идет речь. По сути, величина $\alpha $ играет при этом своего рода роль доверительной вероятности.

Принимая во внимание, что эффективная масса источников гравитационной аномалии определяется достаточно уверенно, возьмем случай $\mu (S*) = \mu (\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} )$. Тогда для $\alpha $ получим следующую оценку: $\alpha \geqslant 1 - {\Delta \mathord{\left/ {\vphantom {\Delta {(2 - \Delta )}}} \right. \kern-0em} {(2 - \Delta )}}$.

Относительно полным результатом математической интерпретации гравитационной аномалии может стать кортеж $\left\langle {{{D}_{1}},{{D}_{2}};S*,\Delta } \right\rangle $.

ПОСТАНОВКА ЗАДАЧИ И ВСПОМОГАТЕЛЬНЫЕ ПОСТРОЕНИЯ

Возьмем традиционные исходные посылки для постановки обратной задачи гравиразведки рудного типа, когда аномалия гравитационного поля $\Delta g$ обусловлена массами, распределенными с постоянными эффективными плотностями ${{\delta }_{k}} > 0$ по оцениваемым связным парциальным носителям ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}}$, $k = 1,2,...,m$, составляющим модель многосвязного возмущающего объекта. Надо подчеркнуть, что такие предпосылки не вполне адекватны реальности, в которой эффективная плотность рудных тел не является строго постоянной. Но такая постановка прижилась в геофизике, и, к тому же, при удачной параметризации модели источников поля эту неадекватность можно снять введением в рассмотрение средних эффективных плотностей ${{\bar {\delta }}_{k}}$ и близких к ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}}$ опорных модельных носителей $S_{k}^{0}$ таких, что невязкой измеренного поля и поля распределения ${{\left\langle {S_{k}^{0},{{{\bar {\delta }}}_{k}}} \right\rangle }_{k}}$ можно пренебречь. Именно опорные носители $S_{k}^{0}$, ассоциируемые с неизвестными истинными телами ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}}$, всегда принимаются по умолчанию за искомые точные решения обратной задачи.

Помимо числа $m$ локальных тел ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}}$ и их эффективных плотностей ${{\delta }_{k}}$, будем считать известными: 1) ограниченные области пространства $D_{k}^{ + }$, содержащие носители ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}}$; 2) области $D_{k}^{ - }$, являющиеся фрагментами носителей ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}}$; 3) оценку $\varepsilon $ максимально возможной нормы помехи в измерениях гравитационного поля; 4) численную характеристику гладкости границы каждого парциального носителя, которую удобно задать в терминах структурных элементов модельного класса $M$ источников поля. Возможны и другие ограничения, используемый нами метод поиска приемлемых вариантов интерпретации их допускает. Задача состоит в построении неулучшаемых (или близких к ним) двухсторонних (множественных) оценок для парциальных носителей ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}}$:

Априорными приближениями к искомым областям $D_{1}^{{(k)}}$ и $D_{2}^{{(k)}}$ являются области $D_{k}^{ + }$ и $D_{k}^{ - }$.

Принимая во внимание базовое представление для областей $D_{1}^{{(k)}}$ и $D_{2}^{{(k)}}$ как объединение и пересечение соответствующих парциальных носителей из всех допустимых модельных носителей, логично предположить, что для построения приемлемых аппроксимаций этих областей понадобится достаточно широкое подмножество допустимых решений обратной задачи в выбранном модельном классе $M$. Это обязывает алгоритм построения таких решений быть полностью формализованным. Таким образом возникают три взаимосвязанные подзадачи: 1) выбор наиболее подходящего модельного класса $M$ источников поля; 2) выбор метода поиска отдельных допустимых решений в классе $M$; 3) выбор алгоритма построения оценок $D_{1}^{{(k)}}$ и $D_{2}^{{(k)}}$, учитывающего особенности класса $M$ и метода отыскания допустимых решений в этом классе.

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

Воспользуемся следующими терминологией и обозначениями. Систему $V$, элементами которой являются связные замкнутые области ${{V}_{j}}$, назовем замощением ограниченной области $D \subset {{{\mathbf{R}}}^{2}}$, если: объединение всех элементов ${{V}_{j}}$ системы содержит область $D$, тогда как общими у двух элементов могут быть разве что их граничные точки; известен алгоритм, позволяющий для каждого ${{V}_{n}} \in V$ установить в заданной системе координат геометрическое место точек $X \in {{V}_{n}}$ и найти номера $j$ всех элементов замощения ${{V}_{j}}$, граничащих с ${{V}_{n}}$. Замощение является регулярным, если все его элементы ${{V}_{j}}$ конгруэнтны некоторому ${{V}_{0}}$. Аналогично определим замощение в R³. Объединение $\Omega $ любого числа $n$ элементов заданного замощения назовем конфигурацией, множество $Ker(\Omega )$ всех элементов замощения ${{V}_{j}} \subset \Omega $ – ядром мощности $n$, а множество $O(\Omega )$ всех элементов ${{V}_{j}} \notin Ker(\Omega )$, граничащих с элементами ядра $Ker(\Omega )$ $ - $ оболочкой конфигурации $\Omega $. В качестве интегральной оценки гладкости конфигурации можно взять отношение мощности ее оболочки к мощности ядра. Модельный класс $M$ носителей, заданный на некотором замощении $V$, определим как множество всевозможных объединений:

(назовем их конфигурациями), составленных из элементов этого замощения.

Перейдем ко второй подзадаче – выбору метода построения отдельных допустимых решений обратной задачи в классе конфигураций (5). Попытки использовать для этой цели методы целочисленного программирования, когда с элементами ${{V}_{j}}$ ассоциируются оптимизируемые бинарные переменные, предпринимались еще в 60-е годы. Помимо проблем вычислительного характера – необходимость решения условно-экстремальных задач большой размерности на маломощных ЭВМ того времени – возникали проблемы принципиального характера, связанные с формализацией требования связности носителя и пофрагментной гладкости его границы.

Подход к построению приближенного решения обратной задачи как к условно-экстремальной задаче, решаемой одним из классических методов минимизации, не единственно возможный. В геофизике (в частности, в геоэлектрике) весьма эффективными зарекомендовали себя алгоритмы интерпретации, основанные на базовой идее метода конечных элементов [Галлагер, 1984] и активно использующие особенности предметной области, в нашем случае – свойства потенциальных полей. В обратных задачах гравиразведки большинство известных алгоритмов такого рода (к примеру, [Zidarov, Zhelev, 1970; Булах, Корчагин, 1978; Rene, 1986; Camacho и др., 2000]) идейно примыкают к монтажному методу Страхова–Лапиной [Страхов, Лапина, 1976], предназначавшемуся его авторами для решения задачи оценивания локального тела $S$ известной плотности $\delta > 0$, удовлетворяющего включению ${{D}^{ - }} \subset S \subset {{D}^{ + }}$. В дальнейшем, преимущественно в работах авторов настоящей статьи, возможности этого метода были существенно расширены, в том числе на случай многосвязного распределения аномальных масс и различных априорных ограничений. Для полноты изложения приведем вкратце суть метода.

Пусть $V$ – выбранное интерпретатором замощение достаточно обширной области пространства, ${{\Omega }^{ - }}$ и ${{\Omega }^{ + }}$ – связные конфигурационные приближения областей ${{D}^{ - }}$ и ${{D}^{ + }}$. Базовый метод Страхова–Лапиной заключается в том, чтобы, отправляясь от конфигурации ${{\Omega }_{0}} = {{\Omega }^{ - }}$, выстроить конечную цепочку локально-оптимальных приближений ${{\Omega }_{i}}$, ${{\Omega }_{{i - 1}}} \subset {{\Omega }_{i}}$, $i = 1,2,...,n$, такую, чтобы массы эффективной плотности $\delta $, распределенные по области $\Omega {}_{n}$, обеспечили приемлемую невязку наблюденного поля. Одна из собенностей метода состоит в том, что плотность $\delta $ переводится в число свободных, оптимизируемых параметров. На шаге $i \geqslant 1$ выбор наилучшего варианта перехода ${{\Omega }_{{i - 1}}} \to {{\Omega }_{i}}$ осуществляется по следующей схеме: 1) из числа пробных носителей ${{\tilde {\Omega }}_{{i,j}}} = {{\Omega }_{{i - 1}}} \cup {{V}_{j}}$, ${{V}_{j}} \in O({{\Omega }_{{i - 1}}}) \cap Ker({{\Omega }^{ + }})$, претендующих на роль очередного приближения ${{\Omega }_{i}}$, путем прямой проверки отбираются те, что не противоречат априорным ограничениям; они составят множество ${{W}_{j}}$ носителей, допустимых к сравнению; 2) для каждого носителя ${{\tilde {\Omega }}_{{i,j}}} \in {{W}_{j}}$ определяется эффективная плотность $\delta _{j}^{*}$, при которой распределение масс $\left\langle {{{{\tilde {\Omega }}}_{{i,j}}},\delta _{j}^{*}} \right\rangle $ обеспечит минимальное значение $\varepsilon _{j}^{*}$ невязки; 3) в качестве локально-оптимального приближения ${{\Omega }_{i}}$ берется носитель $\Omega _{{i,j(i)}}^{*} = {{\Omega }_{{i - 1}}} \cup {{V}_{{j(i)}}}$, на котором (при минимизирующей невязку эффективной плотности $\delta _{{j(i)}}^{*}$) достигается наименьшее значение невязки $\varepsilon _{{j(i)}}^{*}$. Благодаря условию ${{V}_{j}} \in O({{\Omega }_{{i - 1}}})$, связность очередного приближения обеспечивается автоматически, а прием с выведением плотности в число свободных параметров позволяет добиться достаточно малой невязки уже за несколько первых итераций. Из физических соображений ясно, что при некотором $i = n$ выполнится неравенство $\delta _{{j(n)}}^{*} < \delta $, что и станет критерием завершения итерационного процесса.

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

В случае, когда в постановке обратной задачи фигурируют $m$ связных парциальных носителей ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}}$, $D_{k}^{ - } \subset {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}} \subset D_{k}^{ + }$, $k = 1,2,...,m$, по каждому из которых распределены массы эффективной плотности ${{\delta }_{k}} > 0$, метод Страхова–Лапиной притерпевает незначительные изменения. Как и в случае $m = 1$ вся априорная информация, заданная в терминах пространственных областей, формулируется в терминах конфигураций, определенных на выбранном замощении. В частности, области $D_{k}^{ - }$ и $D_{k}^{ + }$ аппроксимируются связными конфигурациями $\Omega _{k}^{ - }$ и $\Omega _{k}^{ + }$, причем конфигурации $\Omega _{k}^{ - }$ служат начальными приближениями ${{\Omega }_{{k,0}}}$ к парциальным носителям ${{\Omega }_{k}}$ в искомом решении $\Omega $ обратной задачи. В работе [Балк и др., 1994] предложен прием, позволяющий на каждой итерации обеспечить одновременный выход всех $m$ подбираемых плотностей в малые окрестности их истинных значений, не увеличивая при этом размерность оптимизационной задачи. Для этого в число оптимизируемых переводится лишь одна из $m$ плотностей (пусть это ${{\delta }_{1}}$), тогда как значения $\delta _{k}^{*}$, $k \geqslant 2$, плотностей ${{\delta }_{k}}$, которые совместно со значением ${{\delta }_{1}} = \delta _{1}^{*}$ минимизируют невязку, согласуются с последним при помощи линейных зависимостей:

Множество пробных вариантов перехода к очередному $i$-му многосвязному приближению формируется за счет элементов замощения ${{V}_{j}} \in Ker(\Omega _{k}^{ + })$, вошедших во все оболочки $О({{\Omega }_{{k,i - 1}}})$ текущих приближений к парциальным носителям ${{\Omega }_{k}}$. Если при выходе из итерационного процесса достигнутая невязка оказалась приемлемой, то допустимое решение $\Omega $ обратной задачи построено.

Чтобы активизировать метод необходимо задаться замощением $V$, что созвучно проблеме выбора размерности модели, которая обычно выводится за рамки метода решения обратной задачи и считается прерогативой самого интерпретатора. Речь идет, прежде всего, о размерах его элементов ${{V}_{j}}$; если класс $M$ построен на регулярном замощении $V$ – то о мере протоэлемента ${{V}_{0}}$. Необходимо, чтобы в классе $M$ нашелся (опорный) элемент $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\Omega } $ такой, что в случае, если расстояние $d(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\Omega } ,\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} )$ не превосходит заданного $\Delta > 0$, оба носителя – $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} $ и $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\Omega } $ – можно было бы считать практически неразличимыми и принять $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\Omega } $ за истинное решение обратной задачи. Как и (в общем случае, многосвязный) носитель $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} $, его модельный представитель $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\Omega } $ останется неизвестным, но здесь важен лишь факт его существования. Идея конечноэлементной аппроксимации источников поля близка к идее сплайн-аппроксимации [Алберг и др., 1972] – в обоих случаях каждый параметр модели отвечает лишь за конкретный участок (будь то среда или область определения функции). Это существенно облегчает выбор дробности замощения в зависимости от значения $\Delta $.

Теперь о вычислительных затратах. Они не так велики, как могут показаться на первый взгляд, к тому же часть вычислений осуществляется в целочисленной арифметике. Объем вычислений практически не зависит от числа $m$ локальных тел ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}}$ в модели источников и определяется в главном суммарным числом элементов замощения ${{V}_{j}}$, которые составят ядро искомого конфигурационного носителя $\Omega $. Метод с очевидностью допускает широкое распараллеливание. Можно вспомнить закон Амдала, хотя ясно и так, что с ростом числа процессоров $p$ зависимость от $p$ времени, затрачиваемого на каждом шаге на отыскание локально-оптимального приближения, будет близка к обратно пропорциональной. За отсутствием совершенного программного продукта мы не можем непосредственно сослаться на опыт решения задачи на многопроцессорных устройствах. Но если за отправную точку взять факт успешного применения метода при решении сугубо трехмерных обратных задач средней сложности на однопроцессорных компьютерах 30-летней давности [Schäfer, Balk, 1993] и связать это с известным законом Мура об удваивании производительности процессоров каждые два года, то легко подсчитать, что на сегодняшний день, по меньшей мере, для задач средней сложности, возможность построения допустимых решений, исчисляемых сотнями, а скорее даже тысячами, не должна вызвать сомнения.

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

непосредственно использовать базовые представления невозможно – множество $Q$ в классе конфигурационных носителей $M$ хотя конечно, но, все же, слишком велико. Однако помимо (8) существуют иные представления областей $D_{1}^{{(k)}}$ и $D_{2}^{{(k)}}$. Из общих соображений ясно, что для каждой пары ($l,k$), $l = 1,2$, $k = 1,2,...,m$ найдутся (причем в большом количестве!) достаточно узкие подмножества $Q(l,k) \subset Q$ решений ${{\Omega }_{r}} = ({{\Omega }_{{r,1}}},{{\Omega }_{{r,2}}},...,{{\Omega }_{{r,m}}})$ таких, что:

Мощность некоторых из этих подмножеств не будет превосходить мощность ядер парциальных носителей ${{\Omega }_{{r,k}}}$. Чтобы построить одно из таких подмножеств необходимо, чтобы каждый очередной из найденных допустимых носителей был результативным, то есть уточняющим текущее приближение к искомой области $D_{l}^{{(k)}}$. Для этого, в свою очередь, необходимо, чтобы каждое допустимое решение обратной задачи было построено при дополнительном ограничении, учитывающем свойства ранее найденных решений.

ПОСТРОЕНИЕ ВЕРХНИХ ОЦЕНОК

Алгоритм построения верхних оценок $D_{1}^{{(k)}}$, $k = 1,2,...,m$, за начальные приближения $D_{{1,0}}^{{(k)}}$ к которым можно взять аппроксимации $\Omega _{k}^{ - }$ априори заданных областей $D_{k}^{ - }$, представляет последовательность самостоятельных, но, вместе с тем, тесно взаимосвязанных итерационных процессов ${{A}_{k}}$, $k = 1,2,...,m$. Задача, которая ставится перед каждым из процессов ${{A}_{k}}$, заключается в обнаружении любого одного из подмножеств $Q(1,k)$. На практике речь идет о формировании подсистемы $\Psi (1,k)$ результативных решений обратной задачи (их число $n(1,k)$ устанавливается апостериори), которых хватит для построения области $D_{1}^{{(k)}}$. С одной оговоркой: хватит с точностью до последствий, связанных с нерешенной проблемой глобального минимума многоэкстремальных функций. В различных приложениях эту проблему принято считать “неизбежным злом”, помнить о ее существовании, но каждый раз не заострять на этом внимание. Мы так и поступим, и вернемся к проблеме надежности в последнем разделе статьи.

Для удобства используется сквозная нумерация итераций, единая для всего алгоритма. На каждой итерации $j$ алгоритма, пришедшейся на работу процесса ${{A}_{k}}$, $k = 1,2,...,m$, осуществляется (одноразовая) попытка найти результативное допустимое решение обратной задачи $\Omega = ({{\Omega }_{1}},{{\Omega }_{2}},...,{{\Omega }_{m}})$, в котором парциальный носитель ${{\Omega }_{k}}$ позволил бы уточнить предшествующее приближение $D_{{1,j - 1}}^{{(k)}}$ к искомой оценке $D_{1}^{{(k)}}$. Для этого достаточно, чтобы ядро носителя ${{\Omega }_{k}}$ содержало хотя бы один элемент замощения ${{V}_{{t(j)}}}$, не вошедший в ядро конфигурации $D_{{1,j - 1}}^{{(k)}}$. Это требование будет выполнено, если выбранный элемент ${{V}_{{t(j)}}}$ включить в ядро нулевого приближения к парциальному носителю ${{\Omega }_{k}}$ без права изымать его из всех последующих приближений. Попытка построения такого решения может оказаться неудачной. Это произойдет тогда, когда допустимых носителей $\Omega $, в которых ядро $Ker({{\Omega }_{k}})$ содержало бы элемент ${{V}_{{t(j)}}}$, попросту нет, либо метод “споткнулся” о проблему глобального экстремума. Однако элемент ${{V}_{{t(j)}}}$ может по воле случая оказаться в ядрах соответствующих парциальных носителей в последующих найденных допустимых решениях, в том числе и тех, что еще будут построены с помощью процесса ${{A}_{k}}$. Поэтому, если на шаге допустимое решение не удалось построить, элемент ${{V}_{{t(j)}}}$ будет лишь временно идентифицирован как не принадлежащий ядру $Ker(D_{1}^{{(k)}})$; окончательно статус этого элемента прояснится по завершению построения оценок (7).

Через $\Omega _{{r(k,i),s}}^{*}$ обозначим $s$-ый парциальный носитель, входящий в $i$-ое по счету результативное решение $\Omega _{{r(k,i)}}^{*}$, построенное с помощью итерационного процесса ${{A}_{k}}$ с учетом дополнительного требования относительно выбранного элемента замощения ${{V}_{{t(j)}}}$. Предположим, что очередной $j$-й шаг алгоритма, пришедшийся на итерационный процесс ${{A}_{k}}$, оказался результативным, причем за предыдущие шаги с помощью этого процесса уже удалось построить $i - 1$ результативных решений. В таком случае найденное решение принимается за очередной, $i$-ый элемент $\Omega _{{r(k,i)}}^{*}$ системы $\Psi (1,k)$ и участвует в построении оценки $D_{1}^{{(k)}}$. Другие $m - 1$ парциальных носителей ${{\Omega }_{s}}$, $s \ne k$, – и в этом суть взаимных связей рассматриваемых итерационных процессов – могут попутно содействовать построению оценок $D_{1}^{{(s)}}$, $s \ne k$:

(предсказать заранее, для каких $s$ окажется, что $Ker(D_{{1,j}}^{{(s)}}) \ne Ker(D_{{1,j - 1}}^{{(s)}})$, невозможно). Если $j$-й шаг не привел к допустимому решению, полагаем $D_{{1,j}}^{{(s)}} = D_{{1,j + 1}}^{{(s)}}$, $s = 1,2,...,m$. Приближения $D_{{1,j}}^{{(s)}}$, $s = 1,2,...,m$, полученные к моменту выхода алгоритма из последнего, $m$-го итерационного процесса ${{A}_{m}}$, принимаются за искомые оценки $D_{1}^{{(s)}}$ с правом подкорректировать их (при необходимости) по результатам работы алгоритма построения нижних оценок $D_{2}^{{(s)}}$.

Открытыми остаются вопросы: как выбирать элемент ${{V}_{{t(j)}}}$ из возможных; каков критерий завершения итерационного процесса; как обеспечить связность носителя ${{\Omega }_{k}}$, если нулевое приближение $\Omega _{k}^{ - } \cup {{V}_{{t(j)}}}$ не является связным.

Начнем с выбора элемента ${{V}_{{t(j)}}}$. Связность области $D_{1}^{{(k)}}$ как объединения конечного числа связных парциальных носителей, имеющих общий связный фрагмент $\Omega _{k}^{ - }$, позволяет сократить число элементов замощения, подлежащих идентификации на предмет принадлежности области $D_{1}^{{(k)}}$. Элементы замощения ${{V}_{{t(j)}}}$, возможно принадлежащие области $D_{1}^{{(k)}}$, но не вошедшие в ядро текущего приближения $D_{{1,j - 1}}^{{(k)}}$, следует искать среди тех, что принадлежат множеству $O(D_{{1,j - 1}}^{{(k)}}) \cap Ker(\Omega _{k}^{ + })$. При некотором $j$ обязательно случится так, что в этом множестве останется лишь один элемент ${{V}_{{t(j)}}}$, который в ходе предшествующих итераций не был опробован на предмет принадлежности ядру $k$-го парциального носителя, причем попытка построить такое решение на итерации $j$ также не увенчалась успехом. Это станет критерием завершения процесса ${{A}_{k}}$.

Формально итерационный процесс поиска допустимого решения можно запустить с любого нулевого приближения и не исключено, что он выйдет на связный допустимый носитель. Но можно действовать надежно и изначально несвязное нулевое приближение дополнить до минимального связного. Организуем итерационный процесс, который позволит элементу замощения ${{V}_{{t(j)}}}$ самостоятельно (в автоматическом режиме) обнаружить конфигурацию $\Omega _{k}^{ - }$:

Процесс завершен, как только на некотором шаге $i = i(0)$ выполнилось условие:

Тогда ядро минимальной конфигурации $\tilde {\Omega }$, связывающей элемент замощения ${{V}_{{t(j)}}}$ с конфигурацией $\Omega _{k}^{ - }$, имеет мощность $i(0)$. Его элементы ${{V}_{{l(i)}}}$ найдем при помощи итерационного процесса, выполняющего обратный ход – от $\Omega _{k}^{ - }$ к ${{V}_{{t(j)}}}$. За элемент ${{V}_{{l(1)}}}$ примем любой элемент ядра $Ker({{\tilde {\Omega }}_{{i(0)}}})$ такой, что

Остальные $i(0) - 1$ элемента ${{V}_{{l(i + 1)}}}$ ядра искомой связующей конфигурации суть любые элементы замощения, удовлетворяющие включению:

ПОСТРОЕНИЕ НИЖНИХ ОЦЕНОК

Без учета результатов построения верхних оценок $D_{1}^{{(k)}}$, за начальные приближения $D_{{2,0}}^{{(k)}}$ к искомым оценкам $D_{2}^{{(k)}}$, $k = 1,2,...,m$, можно взять множества $\Omega _{k}^{ + }$. Приближения $D_{{2,j}}^{{(k)}}$ образуют конечную последовательность сужающихся областей $D_{{2,j}}^{{(k)}} \subseteq D_{{2,j - 1}}^{{(k)}}$. Как и в случае оценок $D_{1}^{{(k)}}$ алгоритм построения нижних оценок $D_{2}^{{(k)}}$ представляет собой цепочку взаимосвязанных итерационных процессов ${{B}_{k}}$, каждый из которых отвечает за соответствующую область $D_{2}^{{(k)}}$ и ориентирован на поиск одного из подмножеств $Q(2,k)$. Попутно каждый процесс ${{B}_{k}}$ может выполнить часть работы, возложенной на процессы ${{B}_{s}}$, $s \ne k$. Фактически его работа заключается в формировании подсистемы $\Psi (2,k)$ результативных решений $\Omega _{{r(k,i)}}^{*}$, $i = 1,2,...,n(2,k)$, для нумерации которых оставим те же индексные обозначения $r(i,k)$, что и в случае процессов ${{A}_{k}}$.

Воспользуемся сквозной нумерацией итераций в пределах алгоритма. На каждой итерации $j$, пришедшейся на работу итерационного процесса ${{B}_{k}}$, предпринимается попытка найти допустимое решение $\Omega = ({{\Omega }_{1}},{{\Omega }_{2}},...,{{\Omega }_{m}})$, в котором парциальный носитель ${{\Omega }_{k}}$ не содержал бы какой-то один элемент замощения ${{V}_{{t(j)}}} \in Ker(D_{{2,j - 1}}^{{(k)}}{\kern 1pt} \backslash {\kern 1pt} \Omega _{k}^{ - })$. Предоставим алгоритму возможность самому выбирать этот элемент случайным образом. После того, как элемент ${{V}_{{t(j)}}}$ выбран из числа возможных, следует наложить запрет на его вхождение в ядра всех приближений к парциальному носителю ${{\Omega }_{k}}$. Начальным приближением к искомому решению служит многосвязная конфигурация ${{\Omega }^{ - }} = (\Omega _{1}^{ - },\Omega _{2}^{ - },...,\Omega _{m}^{ - })$. Попытка построения такого носителя может оказаться безрезультативной хотя бы потому, что допустимых носителей $\Omega $, в которых ядро $Ker({{\Omega }_{k}})$ не содержит данный элемент замощения, попросту нет. Предположим, однако, что итерация $j$ оказалась успешной и до того итерационному процессу ${{B}_{k}}$ удалось найти еще $i - 1$ результативных решений. В этом случае обозначим носитель $\Omega $ как $\Omega _{{r(k,i)}}^{*}$ и осуществим пересчет текущих приближений к оценкам $D_{2}^{{(s)}}$:

(по отношению к каким оценкам $D_{2}^{{(s)}}$, помимо $D_{2}^{{(k)}}$, решение $\Omega _{{r(k,i)}}^{*}$ фактически окажется результативным, предвидеть заранее невозможно). Если $j$-й шаг не привел к допустимому решению, полагаем $D_{{2,j}}^{{(s)}} = D_{{2,j - 1}}^{{(s)}}$, $s = 1,2,...,m$.

Если процесс ${{B}_{k}}$ исчерпал все возможности выбора элемента замощения ${{V}_{{t(j)}}}$, он завершается. Приближения $D_{{2,j}}^{{(s)}}$, $s = 1,2,...,m$, полученные к моменту выхода из последнего, $m$-го итерационного процесса ${{B}_{m}}$ принимаются за оценки $D_{2}^{{(s)}}$, которые удалось построить с помощью одних лишь итерационных процессов ${{B}_{s}}$.

Алгоритм получит логическое завершение, если допустимые носители, найденные с помощью процессов ${{B}_{k}}$, $k = 1,2,...,m$, будут суммированы с построенными областями $D_{1}^{{(k)}}$ с целью возможной корректировки последних, точно также, как построенные оценки $D_{2}^{{(k)}}$ могут быть подправлены с помощью допустимых решений, обнаруженных итерационными процессами ${{A}_{s}}$.

СОСТОЯТЕЛЬНОСТЬ АЛГОРИТМА: ВЫЧИСЛИТЕЛЬНАЯ СЛОЖНОСТЬ И ДОСТОВЕРНОСТЬ ПОСТРОЕННЫХ ОЦЕНОК

Если говорить о трудностях, сопутствующих аналитическому построению оценок (1) и алгоритмическому – оценок (2), то в случае первых основные проблемы носят математический характер (особенно когда в постановке задачи присутствует разнообразная априорная информация), тогда как в случае вторых – вычислительный. Исходя из мощности ядер конфигураций $\Omega _{k}^{ + }$ можно дать, пусть и грубые, оценки сверху для $n(1,k)$. Однако сказать, что построение оценок (7) сводится к поиску

допустимых решений обратной задачи, было бы слишком упрощенно. Если бы поиск допустимых решений мы доверили случаю, сведя все к простому их накоплению без учета свойств уже найденных решений, ни о какой достоверности результатов, под которой мы понимаем соответствие оценок (7) фактическим соотношениям между областями $D_{1}^{{(k)}}$, $D_{2}^{{(k)}}$ и истинными парциальными носителями, речи не могло быть. Принципиально важно, что в построении оценок $D_{1}^{{(k)}}$ и $D_{2}^{{(k)}}$ участвуют результативные допустимые решения, найденные с помощью специально структурированного алгоритма.

Основным фактором, определяющим быстродействие алгоритма и достоверность выполненных построений, являются свойства базового элемента алгоритма – монтажного метода поиска отдельного допустимого решения обратной задачи. Можно предложить достаточно простой способ оценки максимального числа операций, требуемых на реализацию метода, причем эта оценка уточняется в ходе выполнения метода. Метод не относится к классу, так называемых, жадных алгоритмов, когда локальная оптимальность на каждом отдельном шаге обеспечивает оптимальность метода в целом. Что касается надежности вывода о существовании допустимых решений при заданных ограничениях, который будет сделан по результатам работы метода, то здесь можно отметить три позитивных момента. Первый из них в том, что перед методом не ставится (как в методе подбора) задача поиска того единственного решения, которое минимизирует невязку; требуется отыскать любое допустимое решение, обладающее не слишком обременительным дополнительным свойством, либо установить, что таковых нет. Второй позитивный момент заключается в том, что метод допускает различные (не выводящие время его реализации за рамки реального) усложнения структуры итерационного шага, при котором удается еще более сузить интервал неопределенности между достигнутым значением невязки и ее глобальным минимумом. Неожиданным союзником – и это третий позитивный момент $ - $ становится приближенная эквивалентность, предоставляющая методу широкие возможности для выбора. Надо сказать, что использование эквивалентности для геофизики не ново и было одним из первых продемонстрировано В.И. Ароновым [Аронов, 1976].

Свой вклад в повышение быстродействия теперь уже самого алгоритма и, в особенности, надежности итогового результата его применения, вносит сама структура алгоритма построения оценок $D_{1}^{{(k)}}$ и $D_{2}^{{(k)}}$. Речь идет как о связях, существующих внутри каждого процесса ${{A}_{k}}$ и ${{B}_{k}}$ (когда результаты итерации $j$ могут быть подкорректированы на более поздних итерациях), так и о попарных связях итерационных процессов ${{A}_{k}}$ и ${{B}_{k}}$. В последнем случае имеются в виду как связи внутри каждой из двух групп итерационных процессов, так и “межвидовые” связи процессов ${{A}_{k}}$ с процессами ${{B}_{k}}$, которые мы на время обсуждения объединим в одну группу. Помощь, которую каждый из $2m$ итерационных процессов может оказать другим $2m - 1$ процессам, есть смысл разделить на две категории. В первом случае речь идет о содействии, которую один итерационный процесс может получить от других до того, как начнет свою работу, и тогда можно говорить о вкладе в экономизацию вычислений. Во втором случае речь идет о содействии, которое итерационный процесс получает от остальных по окончании своей работы, и здесь уже можно говорить о повышении надежности построений за счет “дружественных” процессов.

Обратимся к фактору экономизации. К моменту входа алгоритма построения верхних оценок $D_{1}^{{(k)}}$ в итерационный процесс ${{A}_{k}}$, $k \geqslant 2$, (это произойдет на итерации $j(0) = n(1,1) + n(1,2) + ... + n(1,k - 1)$) предшествующие ему процессы ${{A}_{s}}$, $s < k$, сформируют для области $D_{1}^{{(k)}}$ нулевое приближение

и лишь итерационному процессу ${{A}_{1}}$ придется начинать работу с нулевого приближения $D_{{1,0}}^{{(1)}} = \Omega _{k}^{ - }$. В еще более выгодном положении находится алгоритм построения нижних оценок $D_{2}^{{(k)}}$. Если воспользоваться решениями ${{\Omega }_{r}} \in Q$, найденными в ходе построения верхних оценок $D_{s}^{{(1)}}$, $s = 1,2,...,m$, то к моменту входа в процесс ${{B}_{k}}$ на шаге $j(0) = n(2,1) + n(2,2) + ... + n(2,k - 1)$ все предшествующие процессы сформируют для процессов ${{B}_{k}}$ следующие нулевые приближения где множество в первых скобках есть вклад итерационных процессов ${{A}_{s}}$.

Совместный вклад в построение оценок $D_{1}^{{(k)}}$ и $D_{2}^{{(k)}}$ “дружественных” итерационных процессов ${{A}_{s}}$ и ${{B}_{s}}$, $s \ne k$ может быть настолько высок, что необходимость в выполнении итерационных процессов ${{A}_{k}}$ и ${{B}_{k}}$ может и вовсе отпасть. Это случится, когда в ходе выполнения итерации $j$, пришедшейся на итерационный процесс ${{A}_{s}}$,$s < k$, окажется, что $D_{{1,j}}^{{(k)}} = \Omega _{k}^{ + }$, либо еще до начала итерационного процесса ${{B}_{k}}$ будет зафиксирован исход $D_{{2,j}}^{{(k)}} = \Omega _{k}^{ - }$. Впрочем, такой исход возможен и раньше, когда пересечение соответствующих парциальных носителей в допустимых решениях, построенных в ходе выполнения итерационных процессов ${{A}_{s}}$, совпадет с конфигурацией $\Omega _{k}^{ - }$.

Чтобы ответить на вопрос, вносят ли “дружественные” итерационные процессы определенный вклад не только в экономизацию вычислений, но и в надежность построений, за которые несут ответственность другие процессы, надо попросту построить оценки $D_{1}^{{(k)}}$ и $D_{2}^{{(k)}}$ силами одних лишь итерационных процессов ${{A}_{k}}$ и ${{B}_{k}}$, а затем сравнить их с прежними оценками. Но все же, реальный интерес представляет сопоставление оценок $D_{1}^{{(k)}}$ и $D_{2}^{{(k)}}$, полученных после выхода из обоих алгоритмов, с приближениями $D_{{1,j}}^{{(k)}}$ и $D_{{2,j}}^{{(k)}}$, с которыми алгоритмы выходили из итерационных процессов ${{A}_{k}}$ и ${{B}_{k}}$. Серийные расчеты на моделях средней сложности убеждают, что итерационные процессы ${{A}_{k}}$ и ${{B}_{k}}$ в подавляющем большинстве способны и самостоятельно справиться с поставленными перед ними задачами.

Вопросы достоверности поднимаются разработчиками алгоритмов решения обратных задач гравиразведки весьма неохотно. Объясняется это тем, что при традиционных представлениях результатов интерпретации трудно определить, где изъяны алгоритма, а где следствие ограниченных возможностей собственно геофизического метода, обусловленные приближенной эквивалентностью и недостатком априорной информации. Приемлемое качество приближенных решений обратной задачи на удачно подобранных модельных примерах нельзя считать надежным критерием состоятельности алгоритма. В определенных классах интерпретационных задач достаточно успешным оказался подход к организации численного эксперимента по оцениванию состоятельности метода, когда часть априорных данных переводится в разряд контрольных, с которыми затем сопоставляются результаты интерпретации, выполненные по оставшимся данным. Пример успешного применения такого подхода при аналитической аппроксимации геопотенциальных полей демонстрируется, в частности, в работах [Ягола и др., 2014; Степанова и др., 2017]. Если такой подход перенести на задачи количественной интерпретации гравиметрических данных, то речь может идти об изъятии из постановки обратной задачи не только части измерений поля, но и отдельных ограничений на параметры модели источников поля, чтобы затем убедиться, удовлетворяет или нет этим ограничениям приближенное решение по неполным данным. При положительном исходе можно допустить, что объем априорной информации достаточен для того, чтобы обеспечить приемлемую близость приближенного и истинного решений. В случае отрицательного исхода вопрос останется открытым. Предложенный в статье алгоритм имеет одно несомненное преимущество – контроль за надежностью процессов построения подмножеств результативных решений обратной задачи здесь не связан с вынужденным изъятием части априорной информации из интерпретируемых данных. “Дружественные” итерационные процессы одновременно выполняют функцию контроля за надежностью алгоритма и, при необходимости, позволяют эту надежность повысить, внося коррективы в основные построения.

Не следует ожидать гарантированного выполнения на практике всего, что обещано в теории, принимая во внимание проблему глобального минимума и, в особенности, неконтролируемую неадекватность модельных посылок. Однако для приложений не столь уж важно, достигнуто или нет полное совпадение построенных областей $D_{1}^{{(k)}}$ и $D_{2}^{{(k)}}$ с объединением и пересечением всех допустимых носителей; мы это все равно не сможем проверить. На практике вполне достаточно, если найденные оценки (7) обеспечили двухсторонние включения для истинных парциальных носителей. Да и здесь, абсолютная безукоризненность построений необязательна, достаточно, чтобы области $D_{1}^{{(k)}}$ и $D_{2}^{{(k)}}$ помогли решить поставленную геологическую задачу. Пародоксальность ситуации состоит в том, что чем сложнее модель источников поля и заметнее проявление приближенной эквивалентности, тем выше у алгоритма шансы обеспечить включения (7), поскольку на каждом шаге ему предоставляется более широкий выбор при поиске очередного результативного решения обратной задачи.

МОДЕЛЬНЫЙ И ПРАКТИЧЕСКИЙ ПРИМЕРЫ

Пример 1. Если, следуя К. Шеннону, понимать под информацией снятую неопределенность об объекте изучения, то мозаичное замощение изучаемой части пространства, порождаемое особенностями взаиморасположения областей $D_{1}^{{(k)}}$, $D_{2}^{{(k)}}$, $D_{k}^{ + }$ и $D_{k}^{ - }$, можно проинтерпретировать следующим образом. 1). Фрагменты $D_{k}^{ + }{\kern 1pt} \backslash {\kern 1pt} D_{1}^{{(k)}}$ не имеют пересечений с носителями ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}}$. 2). Помимо априори известных фрагментов $D_{k}^{ - }$ носителей ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}}$ выявлены новые их фрагменты – области $D_{2}^{{(k)}}{\kern 1pt} \backslash {\kern 1pt} D_{k}^{ - }$; совпадение областей $D_{2}^{{(k)}}$ и $D_{k}^{ - }$ означает, что совокупных данных недостаточно, чтобы обнаружить новые точки носителя ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}}$. 3). Области $D_{1}^{{(k)}}{\kern 1pt} \backslash {\kern 1pt} D_{2}^{{(k)}}$ остаются областями неопределенности, относительно каждой точки каждой из которых нельзя сделать однозначный вывод о (не)принадлежности носителю ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}}$; их меры являются мерой проявления приближенной эквивалентности. 4). Если две или более областей $D_{1}^{{(k)}}$, $k \in K$, имеют непустое пересечение, то его логично назвать зоной повышенного проявления приближенной эквивалентности, каждая точка которой может либо оказаться точкой одного из носителей ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}}$, $k \in K$, либо, вообще, не принадлежать ни одному из них.

Рассмотрим модельную задачу картирования интрузивных тел: диоритового лакколита ${{S}_{1}}$ эффективной плотности 0.15 г/см³ и смещенного по разлому габбрового хонолита, представленного телами ${{S}_{2}}$ и ${{S}_{3}}$ эффективной плотности 0.25 и 0.3 г/см³ соответственно. “Наблюденные” значения (рис. 1а) силы тяжести $\Delta g$ содержат нормально распределенную помеху со среднеквадратическим значением 0.15 мГал, принятым за максимально приемлемое значение нормы невязки “наблюденного” и модельного полей. Помимо обязательной информации заданы некоторые ограничения на максимальную протяженность тел по горизонтали и вертикали (их числовые характеристики сейчас можно не приводить). Чтобы не перегружать рисунок мы также не приводим области $\Omega _{k}^{ + }$ (они были взяты настолько широкими, что не повлияли на построенные оценки), и области $\Omega _{k}^{ - }$ (в качестве последних взяты отдельные элементы замощения).

На рис. 1б приведены стилизованные (со сглаженными границами) области $D_{1}^{{(k)}}$, $k = 1,2,3$, которые можно позиционировать как области, содержащие тела ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{k}}$ (малыми фрагментами интрузий, вышедших за границы области $D_{1}^{{(1)}}$ и $D_{2}^{{(2)}}$, можно пренебречь). Повышенное проявление приближенной эквивалентности не наблюдается. Выявлены две сравнительно небольшие области $D_{2}^{{(1)}}$ и $D_{2}^{{(2)}}$, принадлежащие телам ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{1}}$ и ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{2}}$. Установить какой-либо новый фрагмент тела ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{S} }_{3}}$ не удалось – априорных данных оказалось недостаточно.

Суть множественных оценок проявляется и в их возможности сопоставить по информативности два различных объема априорных данных. На рис. 1в представлены области $D_{1}^{{(k)}}$ и $D_{2}^{{(k)}}$ для случая, когда точных данных об эффективных плотностях нет (известно лишь, что 0.1 $ \leqslant {{\delta }_{1}} \leqslant $ 0.3 и 0.2 $ \leqslant {{\delta }_{2}},{{\delta }_{3}} \leqslant $ 0.3), но заданы более жесткие, чем прежде, ограничения на протяженность тел. Декартово произведение интервалов возможных значений эффективных плотностей покрывалось сеткой, для каждого узла которой определялось “свое” множество допустимых решений обратной задачи. Построенные для всех комбинаций допустимых плотностей ${{\delta }_{k}}$ решения собирались в одно множество, и на основании этого множества строились оценки $D_{1}^{{(k)}}$ и $D_{2}^{{(k)}}$. Как следует из сопоставления рисунков 1б и 1в, объемы извлеченной информации о носителях ${{S}_{k}}$ оказались примерно одинаковы – потери априорной информации об эффективных плотностях компенсируются более жесткими ограничениями на геометрические параметры оцениваемых тел. Правда, здесь уже наблюдается пересечение областей $D_{1}^{{(2)}}$ и $D_{1}^{{(3)}}$, образующее зону повышенного проявления приближенной эквивалентности.

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

Рассмотрим пример, в котором с помощью двухсторонних включений типа (2) надлежит оценить местоположение рудоносной интрузии габбро-долеритов по данным крупномасштабной гравиметрической съемки, выполненной над месторождением медно-никелево-платиновых руд Норильск-1 (рис. 2). Пример интересен тем, что геометрия интрузии довольно хорошо восстанавливается по данным бурения, которые мы не внесли в априорную информацию. Мы воспользуемся ими для контроля. Проследим реакцию оценок на те или иные проявления неадекватности модельных представлений: будет она хаотична или до проведения вычислений можно было из общих соображений предвидеть, как неадекватность скажется на мере областей ${{D}_{1}}$ и ${{D}_{2}}$ и какие подвижки в пространстве с ними произойдут. С этой целью в интерпретационной модели намеренно проигнорирована возможность существования ореола уплотнения вмещающих эффузивных пород, что весьма характерно для надинтрузивных зон месторождений Норильского района. Также взят более жесткий, чем следовало, показатель гладкости границы интрузии.

Приняты следующие допущения. Аномалия $\Delta g$ в основном обусловлена рудоносной интрузией базит-гипербазитового состава, обладающей эффективной плотностью 0.2 г/см³ по отношению к вмещающим породам туфолавой толщи. Предполагаемый уровень помех составляет 0.15 мГал, что примерно отвечает точности гравиметрической съемки. В качестве области ${{\Omega }^{ - }}$, заведомо принадлежащей интрузии, взят отдельный элемент замощения.

Выход за пределы области ${{D}_{1}}$ маломощного фрагмента рудоносной интрузии в левой части геологического разреза объясняется чрезмерно жесткими ограничениями на гладкость границы носителя аномальных масс. Следствием этих ограничений является и несколько завышенная мера области ${{D}_{2}}$. Ее смещение в верхнюю часть разреза относительно фактического положения рудоносной интрузии можно объяснить неучтенным в интерпретационной модели ореолом уплотнения вмещающих эффузивных пород.

И последнее. Независимо от состава и мощности множества допустимых решений, участвующих в построении оценок $D_{1}^{{(k)}}$ и $D_{2}^{{(k)}}$, последние можно трактовать следующим образом: при адекватности модельных представлений любая точка области $D_{1}^{{(k)}}$ может оказаться точкой парциального носителя ${{S}_{k}}$, а фрагменты пространства, которые гарантированно можно было бы отнести к фрагментам этих носителей, малы настолько, что лежат внутри областей $D_{2}^{{(k)}}$. И это уже строгий, геофизически содержательный результат.