CC BY
71
ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА
2010 Геология Вып. 1 (9)
УДК 550.831
Эффективность применения многопроцессорных вычислительных систем с целью оценки достоверности решения обратной задачи гравиметрии
П.И. Балка, А.Г. Деменевь, А.С. Долгальь, О.В. Леденцовь, А.В. Мичурин0
аГермания, Берлин
ьПермский государственный университет, 614990, г. Пермь, ул. Букирева, 15 E-mail: a-demenev@psu.ru, dolgal@mi-perm.ru Торный институт УрО РАН, 614007, г. Пермь, ул. Сибирская, 78 А E-mail: ami.perm@mail.ru
(Статья поступила в редакцию 6 октября 2010 г.)
Рассматривается метод оценки достоверности решения нелинейной обратной задачи гравиметрии (ОЗГ) на основе гарантированного подхода, позволяющий выделять области геологического пространства, заведомо содержащие аномалиеобразующие объекты. Концепция гарантированного подхода предусматривает построение некоторого репрезентативного конечного подмножества допустимых решений ОЗГ. Отмечается необходимость применения высокопроизводительных вычислений для решения гравиметрических задач в 2D и 3Б-вариантах. Отмечаются перспективы создания параллельных технологий решения подобных 3D-задач на базе современных многопроцессорных компьютеров.
Ключевые слова: гравиметрия, обратная задача; оценка достоверности; монтажный метод; гарантированный подход; распараллеливание, высокопроизводительные вычисления, многопроцессорные вычислительные системы.
В самых общих чертах рассматриваемый итерационный процесс становится трехуров-метод оценки достоверности решения нели- невым, в самом внешнем цикле которого нейной обратной задачи гравиметрии (ОЗГ) осуществляется переход к все более размы-представляет собой двухуровневый итераци- тым априорным представлениям о параметрах онный процесс, во внутреннем цикле которо- интерпретационной модели, вносящих ука-го продуцируются эффективные (с точки зре- занную неопределенность. Платой за услож-ния решения основной проблемы) допусти- нение вычислительного процесса является мые решения ОЗГ, а во внешнем - осуществ- возможность упорядочения отдельных подоб-ляется корректировка текущих приближений ластей из D1 по вероятности обнаружения в D1 и D2 к искомым областям D1 и D2, пер- их пределах источников гравитационной ано-вая из которых содержит аномальные массы, малии.
заключенные в неизвестном объеме ST, а Центральное место в техн°л°гиях п°-другая целиком заполнена фрагментами ано- строения пары < Dl, D2 > занимает собст-
мальных масс: D2 с ST с D1. Если априор- венно алгоршм построения отдетмьк допустимых решений. От него требуется не только ная информация содержит (как это всегда и ^ ^ J
высокая эффективность по быстродействию,
бывает на практике) неопределенность, то
© Балк П.И., Деменев А.В., Долгаль А.С., Леденцов О.В., Мичурин А.В., 2010
50
но и способность вырабатывать «особые» допустимые решения, небольшое число которых могло бы составить некое репрезентативное подмножество, коллективно обладающее теми же свойствами, что и все множество допустимых решений (являющееся, вообще говоря, бесконечным). В плане обеспечения приемлемого времени построения отдельно взятого допустимого решения весьма желательно, чтобы алгоритм допускал распараллеливание процесса вычислений, обеспечивающее возможность использования многопроцессорной вычислительной техники. Говоря о проблеме целенаправленного поиска сравнительно небольшого числа допустимых решений ОЗГ, обладающих свойством репрезентативности, можно утверждать, что решить ее можно, лишь опираясь на идею аппроксимации.
Сама по себе эта идея не нова. Неявно она присутствует во всех известных алгоритмах решения ОЗГ, где по умолчанию всегда принимается, что природное распределение масс можно с приемлемой точностью описать элементом модельного класса с небольшим числом свободных параметров. Как методообразующая идея аппроксимации наиболее ярко выражена в технологиях, использующих дискретное (сеточное) описание плотностной среды, когда объем (носитель) аномалиеобразующих масс представляется в виде объединения некоторого числа элементарных геометрических фигур (к примеру, призм). Подобные технологии построения решений ОЗГ представлены давно уже сданными в архив палетками из эпохи ручного счета. Пик наиболее интенсивных работ по созданию компьютерных технологий решения ОЗГ, базирующихся на использовании сеточных моделей (ему - отдадим должное - предшествовали работы А.А. Юнькова и Е.Г. Булаха), приходится на конец 60 - начало 70-х гг. (исследования С. В. Шалаева, В. В. Ломтадзе, В.А.Дядюры, Л.Г. Перфильева), и обусловлено это было тем, что многие геофизики-теоретики овладели к тому времени методами линейного программирования, очень удачно вписывающимися в проблематику линейных обратных задач по оценке закона распределения плотности масс в заданном объеме. Попытки использовать методы линейного (в том числе и целочисленного) программирования при решении нелинейной ОЗГ по оценке геометрии области, занятой массами известной (постоянной) плотности, не увенчались успехом. Причиной тому были даже не слабые
возможности вычислительных средств того времени, сколько плохая совместимость сеточных моделей и классических методов минимизации в вопросах учета разнообразной априорной информации.
В 70-х гг. отечественными геофизиками проводились масштабные исследования, смысл которых заключался в поиске наиболее эффективных (среди известных) методов решения условно-экстремальных задач для целей минимизации невязки в процессе подбора допустимых решений нелинейной обратной задачи в раличных модельных классах, отличных от сеточных. Но известные общематематические методы минимизации не привязаны к специфике какого-то определенного класса практических задач, а универсальность, как известно, всегда входит в противоречие с эффективностью. Потому неудивительно, что отыскать подходящий на все случаи алгоритм (или хотя бы четко очертить классы интерпретационных моделей, для которых оптимальным является тот или иной метод минимизации) так и не удалось [1,2]. Зримо обозначилась необходимость разработ -ки в недрах математической геофизики собственных проблемно-ориентированных методов минимизации, жестко привязанных к специфике обратных задач гравиметрии и магнитометрии.
В развитие этой идеи с небольшим интервалом по времени выходят работы А.В. Овча-ренко [з] и В.Н. Страхова [4] (в соавторстве с М.И. Лапиной), в которых предлагался, по существу, новый принцип структурирования итерационного процесса подбора допустимого решения нелинейной ОЗГ в сеточных классах источников поля, обеспечивающий удобный контроль за соблюдением требования связности и односвязности подобранного приближения к неизвестному носителю аномалиеобразующих масс. В.Н. Страхов предложил называть его монтажным. Конструктивно этот принцип был реализован им в алгоритме регулируемой направленной кристаллизации (РНК). В последующем различные аспекты монтажных алгоритмов затрагивались в работах И.И. Корчагина, У. Шеффера, Т.В. Балк, А.С Долгаля. В цикле работ П.И. Балка метод РНК был обобщен, а сам монтажный подход распространен на постановки обратных задач магнитометрии, случай многосвязной аномалиеобразующей среды с различными значениями плотностей парци-
альных источников, на задачи совместной интерпретации различных компонент гравитационного и магнитного полей, на случай смешанной структурно-рудной задачи, наконец - на расширение типов априорной информации, доступных одновременному учету при построении допустимых решений. Справедливости ради надо отметить, что теоретические разработки по рассматриваемой проблематике существенно опережают их практическую апробацию.
Сущность монтажного подхода - а именно он предлагается авторами настоящей статьи в качестве базового элемента технологии оценки достоверности результатов интерпретации гравитационной аномалии - неразрывное единство конечно-элементного способа описания плотностной среды и специального, и что принципиально важно, геофизически-содержательного (монтажного) принципа структурирования итерационного процесса собственно построения приближенного решения. В монтажных технологиях используются следующие основные понятия. Пусть 8
- достаточно широкая область нижнего полупространства, заведомо содержащая все источники аномалии, сосредоточенные в объеме
(носителе) 8Т. Конечно-элементным (сеточным) замощением области 8 назовем представление последней в виде объединения конечного числа «элементарных» замкнутых подобластей (геометрических фигур) о а,
любые две из которых имеют общие разве что граничные точки. В двумерных постановках ОЗГ это могут быть квадраты или шестиугольники, в трехмерной - параллелепипеды; возможна и более сложная структура элемента, учитывающая часть априорной информации (к примеру, предполагаемые углы наклона аномалиеобразующих тел). Под монтажными классами источников поля понимают множество всевозможных распределений тяготеющих масс, размещенных с заданными плотностями по конфигурационным областям
П = П( А) = , ае А, составленным
из произвольного числа монтажных элементов оа и удовлетворяющим определенным
априорным ограничениям, прежде всего - условиям связности и односвязности. Основные операции в сеточных классах выполняются с использованием понятий ядра Я [о] оболочки о[п] и границы Г [о] конфигурации П:
Я [О] - суть множество элементов Юа е П; О[П] - множество всех элементов замощения Юа , не принадлежащих ядру Я [о] , но граничащих хотя бы с одним из его элементов; Г [п] - множество всех элементов ядра конфигурации П, граничащих с элементами ее оболочки.
В постановке обратной задачи для изолированного тела 8Т известной постоянной плотности оТ > 0 монтажный принцип решения состоит в том, чтобы, отправляясь от
заданной связной конфигурации П0, выстроить конечную последовательность
П0, П1, П2,...,конфигурационных носителей
и ассоциированную с ней последовательность
* *
о1 ,о2,..., плотностей, минимизирующих невязку, таких, что на некотором шаге п поле
*
масс плотности оп, распределенных по конфигурационной области Пп, создавало бы поле, удовлетворяющее наблюденному по невязке, сам конфигурационный носитель Пп удовлетворял бы всей априорной информации
о неизвестном носителе 8Т, а значение плотности оп было бы близко к значению истин-
Т
ной плотности о .
В рамках монтажного принципа можно строить целое иерархическое семейство алгоритмов, различающихся структурой (степенью сложности) итерационного шага. В самом простейшем случае - методе регулируемой направленной кристаллизации (РНК) - любое очередное приближение Пп образуется путем внесения в ядро Я [пп-1 ] одного элемента опа
из О[пп-1 ], обеспечивающего наименьшую среднеквадратическую невязку наблюденного и модельного полей при некотором подобранном значении плотности пробной конфигурации оп > оТ. Распараллеливание допускает уже отдельный шаг < Пп-1, о п_1 > ®
< Пп,Оп > . Критерием достижения завершения итерационного процесса является выполнение условия о п » оТ . Если при этом значение достигнутой невязки подобранного и наблюденного полей укладывается в представления интерпретатора о предполагаемом уровне случайных помех и помех геологического происхождения, то построенная конфи-
гурация Пп принимается за допустимое приближенное решение П рассматриваемой нелинейной ОЗГ.
Суммируя сказанное в адрес монтажных алгоритмов решения ОЗГ, можно утверждать, что монтажный подход снимает проблему неустойчивости приближенных решений обратной задачи в ее классическом понимании, т. к. конечная размерность модели и природные ограничения на ее плотностные и геометрические параметры изначально приводят к компактному множеству возможных решений, а учет определенного объема априорной информации обеспечивает содержательность решения поставленной геологической задачи. Плюс к этому - монтажные алгоритмы чрезвычайно приспособлены к одновременному учету разнообразной априорной информации, характеризующей местоположение, форму и размеры аномалиеобразующих тел; весь вопрос здесь в контроле за сохранением текущим приближением свойств предшествующего приближения к искомому решению, что особо просто осуществляется в методе РНК, где изменение ядер текущих приближений минимальное от итерации к итерации.
Авторы полагают, что изложенные выше соображения достаточно убедительно свидетельствуют в пользу использования монтажного подхода в алгоритмах оценки достоверности решения обратных задач. Теперь такой алгоритм выглядит следующим образом.
Для заданной совокупности априорных ограничений (либо, если некоторые ограничения носят размытый характер, то для фиксированных интерпретатором априорных допущений, укладывающихся в рамки существующей неопределенности) с помощью монтажного алгоритма строится последовательность допустимых решений П*,П*2,..., где каждое последующее решение П*п+1 содержит хотя бы один элемент замощения о а , не вошедший в ядра уже построенных носителей п* , к = 1,2,...,п . Чтобы обеспечить это условие, достаточно итерационный процесс построения допустимого носителя Пп+1 начать с конфигурации (нулевого приближения), содержащей какой-либо элемент оа оболочки объединения уже построенных допустимых конфигурационных носителей П* (назовем этот элемент порождающим), и не исключать
его из ядер последующих приближений, если структура итерационного шага такое принципиально не исключает. Итерационный процесс построения допустимых решений завершается, как только на некотором шаге п +1 не удается более построить допустимое решение при любом выборе порождающего элемента Юа . Тогда семейство конфигурационных носителей п*, п*2,...,п*п считается репрезентативным, объединение построенных допустимых носителей принимается за оценку П* искомой области П, а их пересечение
- за оценку П2* области П2.
Возьмем, наконец, самый общий случай, когда априорная информация несет в себе элемент неопределенности и реализация алгоритма оценки достоверности результатов интерпретации сводится к трехуровневому итерационному процессу. Для наглядности рассмотрим типичную ситуацию, когда уровень (норму) е помех невозможно точно оценить достоверно, известно лишь, что неизвестное значение е лежит в некотором интервале
[еш1п,ешах ]. Возьмем два любых значения предполагаемого уровня помех - е1 и е2 > е1 - и пусть области, гарантированно содержащие (при этих предположениях и прочей априорной информации) есть ) и П22) соответственно. Очевидно, что всегда П^2) ^ . Тогда можно утверждать: веро-
ятность того, что часть аномальных масс заполняет объем П^2), выше вероятности
встречи аномальных масс в объеме .
Ясно, что интервал возможных значений нормы помехи можно по желанию покрыть сколь угодно густой сеткой значений е и построить детальное разбиение области, заведомо содержащей источники поля, по вероятности встречи аномальных масс в каждом из элементов этого разбиения.
На рис. 1 представлены результаты интерпретации аномалии Дg, обусловленной двумя однородными призмами плотности 0,15
3
и 0,1 г/см , при различных предположениях относительно уровня помех в измерениях. Неизвестными являются параметры положения тел в пространстве. Как и следовало ожидать, существует определенный диапазон уровня помех, в рамках которого гарантиро-
ванно удается выявить фрагмент лишь одного (залегающего менее глубоко) возмущающего тела (см. случай £=1,2 мГал).
На рис. 2 приведены результаты решения ОЗГ по данным крупномасштабной гравиметрической съемки, выполненной над месторождением платино-медно-никелевых руд Норильск-1. В основу модели среды положены следующие допущения: аномалия в основном обусловлена рудоносной интрузией базит-гипербазитового состава, ее избыточная плотность (по отношению к вмещающим породам трапповой формации) составляет 0,2 г/см. С помощью монтажного метода при различных центрах кристаллизации (начальных приближениях) было построено около 400 различных вариантов пространственных распределений масс, удовлетворяющих априорным допущениям. Их синтез позволил выделить область П2 (рис. 2, контур 8), с высокой степенью достоверности принадлежащей источнику аномалии при предполагаемом уровне помех £ = ±0,15 мГал.
Оценка достоверности результатов количественной интерпретации в терминах пары множеств П1, П2 может существенно повы-
¿д, мГал
16 т
0 км 10км 20км 30км 40км 50км
Рис.1. Результаты интерпретации гравитационного поля, созданного двумя призмами, в терминах пары множеств < 02 > : 1 - график А^; фрагменты среды 02 , гарантированно принадлежащие
источника м при условии, что среднеквадратическая невязка £ полей не превосходит: 2 - 1,2 мГал, 3 - 0,7 мГал, 4 - 0,4 мГал; 5 - аномалиеобразующие призмы; 6 - границы областей , гарантированно содержащие источники аномалии в предположении, что невязка £ не превосходит 0,4 , 1,5 и 2,5 мГал
сить вероятность вскрытия искомых аномалиеобразующих объектов в заданных интервалах глубин поисковыми и разведочными скважинами, рекомендованными по гравиметрическим данным.
Вышепредставленный подход решения обратной задачи гравиметрии и оценки достоверности результатов количественной интерпретации до недавнего времени был программно реализован с использованием системы объектно-ориентированного программирования Delphi 7.0 для Windows. В этой последовательной программе вычисления осуществлялись на одном вычислительном ядре процессора, что приводило к ощутимым временным затратам. Для построения пары множеств D1, D2 требуются сотни - иногда тысячи - допустимых решений ОЗГ. Время последовательного счета даже в 2D-варианте с использованием одного процессора типа Intel Core составляет при 100-150 точках задания поля и всего одном возмущающем объекте более часа, а в отдельных случаях (при сложных задачах) - и сутки.
Рис.2. Результаты интерпретации поля месторождения Норильск-1: 1 - породы туфовой толщи; 2 - отложения тунгусской серии; 3 - силлы габбро-долеритов; 4 - рудоносная интрузия; 5 - дизъюнктивные нарушения; 6 - локальная составляющая наблюденного поля Д£; 7 -контур области 01, содержащей все источники локальной аномалии Д£; 8 - контур области 02, гарантированно содержащей фрагмент аномалиеобразующего объекта; 9 - буровые скважины
При решении подавляющего большинства практических задач требуется применение высокопроизводительных вычислений даже в 2D, а тем более 3Б-варианте. Решено было использовать распределенные вычисления, которые обычно эффективны для решения вычислительно сложных задач, не требующих интенсивного обмена информацией между параллельными подзадачами. Для реализации этих вычислений были использованы средства распределенной вычислительной инфраструктуры программы «Университетский кластер» («УК»). Официальный сайт программы
- unicluster.ru. Участие в этой программе принимают несколько институтов РАН и многие российские вузы, в том числе и Пермский госуниверситет.
В настоящее время на высокопроизводительных многопроцессорных системах Linux установлен значительно чаще чем, Windows, поэтому важно было обеспечить возможность использования компьютеров с различными операционными системами. Была создана переносимая программная реализация алгоритма оценки достоверности решения ОЗГ в 2D-варианте. Для этой цели использовался кросс-платформенный компилятор FreePascal Compiler и совместимая с ним среда разработки Lazarus.
Parallel Compute: MPI - сервис уровня инфраструктуры программы «УК», который обеспечивает разработку и выполнение MPI программ для вычислительных систем с распределенной памятью. Практически на всех высокопроизводительных вычислительных кластерах установлена коммуникационная библиотека стандарта MPI. Поэтому предложена схема распараллеливания вышеуказанной программы на основе MPI.
Одной из характеристик эффективности параллельного алгоритма является ускорение
- отношение между временем выполнения задания с данными размера m асимптотически оптимальным последовательным алгоритмом и временем выполнения того же задания параллельным алгоритмом на машине с p процессорами. Для априорной оценки ускорения применена сетевая формула Амдала, которая учитывает потери времени на межпроцессорный обмен данными в параллельном приложении:
Sp (m) =---------1---г—,
pW D + / + (l -/)/p
где / = /(n) — доля последовательных операций, D = DalgDtech — коэффициент сетевой деградации. При этом Dalg определяет алгоритмическую составляющую коэффици-
ента деградации, обусловленную свойствами алгоритма, а О1есИ — техническую составляющую, зависящую от соотношения технического быстродействия процессора и аппара -
туры сети. Далее, ^а1% Ж сотт^ ^сотр , где
Ж „
сотт — количество операций передачи
Ж
comp
количество вычисли-
D tech t comm ^ t comp где
данных, а тельных операций.
Котт --- среднее время выполнения одной
операции передачи данных, 1сотр — среднее
время одной вычислительной операции.
Рассмотрена типичная 2Б-задача размером т=19200 монтажных элементов. Оценено соответствующее типичное соотношение последовательной и параллельной части алгоритма / = 2,8 -10 5. На учебном вычисли-
тельном кластере компьютерного центра механико-математического факультета ПГУ оценен соответствующий коэффициент сетевой деградации О = р -1,9 -10-б .
Оценка показывает, что применение высокопроизводительных вычислительных кластеров позволяет кардинально сократить время расчета. В результате распараллеливания программы время счета удается существенно уменьшить (максимально в 340 раз при числе используемых процессоров 1024). Дальнейшее увеличение количества вычислительных узлов (процессоров) лишь замедляет работу программы, т. к. накладные расходы на передачу информации между узлами сети начинают превосходить выигрыш за счет распараллеливания (рис. 3).
Рис.3. Зависимость ускорения от числа процессоров
Развитие гарантированного подхода к интерпретации поля силы тяжести во многом обусловлено большими потребностями геофизической отрасли в оперативной интерпретации больших объемов цифровых данных и тесно связано с возросшими вычислительными возможностями компьютеров. В дальнейшем, по нашему мнению, актуальным в области гарантированного подхода будет проведение исследований по следующим направлениям:
1. Параллельные реализации предложенного способа для 3Б-задач.
2. Разработка алгоритма монтажного метода в модификации регулируемой направленной раскристаллизации, разноуровневой перекристаллизации для 2Б и 3Б-задач.
3. Проектирование и реализация в инфраструктуре программы «Университетский кластер» проблемно-ориентированного вебсервиса для автоматизации интерпретации результатов гравиметрии.
Работа поддержана в рамках проекта «Развитие центра коллективного пользования высокопроизводительными вычислительными ресурсами - НОЦ ПиРВ» Программы разви-
тия национального исследовательского уни-верситета ПГУ.
Работа выполнена при поддержке Уральского отделения РАН (по результатам конкур-
са научных проектов молодых ученых и аспирантов 2010 г.) и РФФИ (грант № 10-05-96023-р_урал_а).
Библиографический список
1. Гольдшмидт В.И. Оптимизация процесса количественной интерпретации данных гравиразведки. М.: Недра, 1984. 184 с.
2. Овчаренко А.В. Подбор сечения двухмерного тела по гравитационному полю// Вопросы нефтяной и рудной геофизики. Алма-Ата, 1975. Вып. 2. С. 71-75.
3. Страхов В.Н., Лапина М.И. Монтажный метод решения обратной задачи гравиметрии// Докл. АН СССР. 1976. Т. 227, №2. С. 344-347.
4. Балк П. И. Использование априорной информации о топологических особенностях источников поля при решении обратной задачи гравиметрии. // Докл. АН СССР. 1989. Т. 345, № 5. С. 1082-1084.
5. Балк П.И. Долгаль А. С., Балк Т.В. Сеточные модели плотностной среды и опыт их применения при прослеживании дифференцированных интрузий по данным гравиразведки // Геология и геофизика. 1993. № 5. С. 127 - 134.
The Efficiency of Multiprocessor Computation Systems Using to Estimate the Reliability of the Solution of Gravimetric Inverse Problem
P.I.Balka, A.G. Demenevb, A.S.Dolgalb, O.V. Ledentsovb, A.V.Michurinc
a Germany, Berlin
b Perm State University, 614990, Perm, Bukirev st., 15 E-mail: : a-demenev@psu.ru , dolgal@mi-perm.ru
c Mining Institute of Ural Branch of RAS, Russian Academy of Sciences, 614007, Perm, Sibirskaya st., 78A. E-mail: ami.perm@mail.ru
We consider a problem of the estimation of the reliability of the solution of nonlinear gravimetric inverse problem (GIP) on basis of guaranteed approach that allows to select areas of the geological space, obviously containing anomaly-creating objects. Guaranteed approach is based on the examination of finite sets of feasible solutions of GIP in the assembly class of field sources. We are noted necessary to apply high performance computing for solving large- scale problems, even in 2D, and even more 3D-version. We are noted prospects for parallel technologies of like3D-task solving at modern multiprocessor computers.
Key words: gravimetry, inverse problem, the reliability estimation, the assembly method, guaranteed approach, parallelization, high performance computing, multiprocessor computing systems.
Рецензент - доктор геолого-минералогических наук Б.А. Спасский
