CC BY
96
ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА
2010
Геология
Вып. 1 (9)
УДК 550.831.017
Алгоритмы истокообразной аппроксимации геопотенциальных полей для данных большой размерности
А.Ф. Шархимуллин, А.С. Долгаль
Пермский государственный университет, 614990, Пермь, ул. Букирева, 15 E-mail: dolgal@mi-perm.ru, art-perm@mail.ru (Статья поступила в редакцию 11 октября 2010 г.)
Рассматривается ряд вопросов, связанных с аналитическими аппроксимациями геопотенциальных полей при большом числе точек наблюдений. Приводятся описания эффективных алгоритмов, примеры модельных и практических интерпретационных построений. Предлагается принципиально новый алгоритм построения истокообразных конструкций, основанный на методе конечных элементов. Намечены направления дальнейших исследований в данной области.
Ключевые слова: истокообразная аппроксимация, трансформация, источники поля, квадродерево, вейвлет-анализ, метод конечных элементов.
Представление геолого-геофизических данных в виде различного рода моделей является необходимым этапом процесса интерпретации. В целом ряде задач, связанных с интерполяцией, трансформациями, фильтрацией численных данных, построение цифровой модели источников поля представляет собой промежуточное звено, а не конечный результат. По мнению академика В.Н. Страхова, методообразующая идея аналитической аппроксимации занимает центральное место в единой теории трансформации геофизических полей [6,7]. Аналитические аппроксимации геопотенциальных полей физически содержательными гармоническими функциями, в качестве которых выступают поля элементарных источников (шаров, вертикальных стержней и т. п.), хорошо зарекомендовали себя на практике. Подобные аппроксимации зачастую называют также «истокообразными аппроксимациями» [2], а саму процедуру их построения - построением аналитических моделей геофизических полей [4]. Следует отметить, что методы истокообразной аппроксимации в основном разрабатывалисьсь применительно к гравитационному и магнитному полям в рамках континуальной постановки.
Использование алгоритмов истокообразной аппроксимации позволяет создавать сеточные эквивалентные распределения петро-физических параметров без привлечения до-
полнительной информации о геологической среде. Один из наиболее простых и широко применяемых методов заключается в построении совокупности эквивалентных источников ниже поверхности наблюдений (горизонтальной плоскости или поверхности рельефа) в узлах регулярной сети точек на глубине 1Ах<Н<2Ах, где Ах - шаг сетки. В качестве возмущающих объектов принимаются тела простой формы (шар, тонкий стержень и т.д.). Решение задачи сводится к решению хорошо обусловленной системы линейных алгебраических уравнений (СЛАУ) вида ОЛ=и, где G - матрица, элементы которой представляют собой истокообразные функции
- поля элементарных источников при единичных коэффициентах условных масс, а и -массив исходных значений поля (рис.1). Количество источников обычно соответствует числу точек задания геопотенциального поля [1].
На основе решения СЛАУ определяются параметры элементарных источников при вычислении поправок за влияние рельефа в гравитационном поле [3]. Единожды построенная аналитическая аппроксимация в дальнейшем может многократно использоваться для вы -числения поправок за удаленные области рельефа на других площадях гравиметрических съемок.
© Шархимуллин А.Ф., Долгаль А.С., 2010
Большое значение в процессе интерпрета -ции реальных геофизических данных имеет фактор размерности решаемой задачи.
Рис.1. Классический метод построения аналитических аппроксимаций
Объем информации в пределах одного участка исследований может составлять десятки и сотни тысяч точек задания поля и более. Поэтому крайне важно минимизировать количество источников в модели, сохранив необходимую точность аппроксимации. Для этих целей и были разработаны аппроксима-ционные конструкции, использующие в качестве вспомогательных инструментов методы квадродерева и вейвлет-анализа.
Основная идея квадродерева заключается в представлении исходной области задания поля в виде ранговых блоков методом квадродерева [4, 5]. Первоначально матрица значений поля делится на небольшое число квадратов первого порядка (рис.2). С целью создания “грубого образа” модели эквивалентные источники размещаются в узлах разреженной сети в центрах ранговых областей на глубине 1Ах1<Н1<2Ах1 (Ах1 - шаг сети первого порядка
- расстояние между центрами квадратов). Уточнение аппроксимационной конструкции на следующем этапе производится путем детализации рангового блока на четыре подобласти - квадраты меньшего размера - и созданием сеточных моделей эквивалентных источников более высокого порядка с меньшим шагом Ах2,...,Ахп и соответственно на меньших глубинах кь...,Нт где п - глубина квадродерева или количество уровней детализации модели. Первым приближением при аппроксимации служит поле эквивалентных источников первого, наиболее глубоко расположенного, уровня - региональная состав -ляющая. Разность наблюденного и модельного полей представляет собой объект аппроксимации на следующем этапе подбора параметров источников. Построение модели завершается, когда достигнута максимальная глубина квадродерева п (шаг сети наиболее высокого порядка Ахп равен расстоянию между узлами исходной матрицы поля) или отклонение модельного поля от измеренного на
каком-либо этапе процесса достигает заданной погрешности аппроксимации е.
источники
источники
Рис.2. Схема расположения эквивалентных источников в центрах ранговых блоков методом квадродерева, модельное гравитационное поле семи сфер и план расположения источников в аппроксимационной конструкции
Таким образом осуществляется адаптация геометрии аппроксимационной конструкции к морфологическим особенностям исходного геофизического поля. Результирующая модель представляет собой совокупность уров-ней с уплотненным распределением источников в плане вблизи ярко выраженных аномальных особенностей поля. В случае, если исходное геопотенциальное поле носит спокойный характер, для его описания требуется сравнительно небольшое количество источников, расположенных в узлах редкой сети.
Алгоритм вейвлет-преобразования предназначен для построения аппроксимаций гравитационного поля совокупностями эквивалентных источников при значительном объеме исходных данных. При этом используется эффективный алгоритм «истокообразной аппроксимации», в котором построение адаптированной к полю аппроксимационной конструкции осуществляется на основе предварительно выполненного вейвлет-преобразования наблюденного поля с базисными функциями Хаара (Долгаль и Симанов, 2008).
Рис.3. Сеточные распределения эквивалентных источников геопотенциальных полей с использованием быстрого вейвлет-преобразования
Любой сигнал при проведении быстрого вейвлет - преобразования может быть представлен в виде некоторого «грубого образа» и уточняющих его разномасштабных «уточняющих особенностей» (рис.3), что осуществляется его разложением с использованием ортонормированных масштабируемых и перемещаемых в пространстве функций с компактным носителем вейвлетов. Существуют устойчивые численные алгоритмы вычисления вейвлет-коэффициентов при таком разложении, а также обеспечивается возмож-ность приближенного восстановления сигнала с априори заданной точностью. Поскольку
вейвлет-коэффициенты существенно отличаются от нуля только вблизи сингулярностей (т.е. вейвлет-ряды «нормальных» функций допускают сильное «разряжение»), удается отбросить мелкие детали и выделить наиболее существенные особенности анализируемого сигнала, что позволяет существенно сжать объем первоначальной информации.
Авторами разработан принципиально иной подход к численной процедуре истокообразной аппроксимации гравитационного поля, базирующийся на методе конечных элементов (МКЭ). Ранее метод конечных элементов нашел свое применение при решении задач количественной интерпретации гравитационных и магнитных аномалий (монтажные алгоритмы решения нелинейных обратных задач) [6].
Основной особенностью предлагаемого авторами метода построения аналитических аппроксимаций гравитационного поля является его радикальная декомпозиция - проблема сводится к сумме простейших одномерных задач линейной минимизации. Его концептуальной основой служит идея метода конечных элементов, ранее нашедшего применение при решении задач количественной интерпретации гравитационных и магнитных аномалий (монтажные алгоритмы решения нелинейных обратных задач) [1].
Реализация метода начинается с выбора некоторого (в общем случае нерегулярного)
замощения Т = {0( 1)}] части S геологического пространства, заведомо содержащего источники рассматриваемого поля, объемными геометрическими фигурами 0(1), для которых (при постоянной плотности) существуют достаточно эффективные алгоритмы решения прямой задачи (в наиболее благоприятном случае - аналитические формулы). В частном и наиболее простом случае - это замощение (вообще говоря, разновеликими) параллелепипедами. Объединение
п= и п(,) , |/| = т, I <^{1,2,...,п} элементов некоторого подмножества Т0 ^ Т , взятых с различными (но постоянными) «оптимальными» плотностями о(г) = &°р‘(1), и будет выступать в качестве «фиктивного» источника, полем которого аппроксимируется наблюденное аномальное гравитационное поле и(X) . Так, построенное множество О принято называть конфигурацией, а само под-
множество Т0 элементов замощения, объединенных в множестве О - ее ядром я[о]. Выбранный способ представления аппроксимирующего (по полю!) источника предопределяет необходимость минимизации числа п элементов замощения О(/), и в конечном итоге - мощности т ядра конфигурационного носителя О (по сути, ставится вопрос об экономизации вычислительного процесса). Формализованный критерий выбора минимального числа п, обеспечивающего разумную точность аппроксимации, найти пока не удалось. Ясно, однако, что значение п должно быть согласовано, прежде всего, с масштабом работ - сетью измерений гравитационного поля. Однако уже сейчас, исходя из эвристических соображений (численные расчеты это также подтвердили), можно сказать, что суммарные (по всей сети измерений) гравитационные влияния элементов замощения по возможности должны иметь один и тот же порядок. В качестве ориентира при выборе размеров элементов замощения О(/) может выступить, в частности, условие
N
£и 2(О(0; X ) = се2, (1)
,=1
где N - число точек, в которых задано аппроксимируемое аномальное поле, и (0(/); X,) - гравитационное влияние тела
О(/) единичной плотности в точке X, (х,, у,, Zj) произвольного рельефа, £ -требуемая точность приближения аномального поля в метрике Ь2, с - некоторая константа (предположительно с » 1). Существует несколько способов упорядочения процедуры построения подобных замощений. Можно, в частности, задавшись последовательностью
z < z < < z V
1 2 Ь, рассечь область ° парал-
лельными
процесс, на внешнем шаге которого осуществляется наращивание конфигурационного носителя на один элемент из 7 , а на внутреннем
- собственно выбор «оптимального» элемента замощения (из тех, что не вошли в ядро конфигурации, построенной на предыдущих итерациях внешнего цикла), обеспечивающей минимум невязки.
В развернутом виде это выглядит следующим образом. Пусть
г-1
О-1 = и по, )
к=1
(1(г -1) = {h,/2,...Л-l}^{1,2,...,п}) (2)
- приближение к искомой конфигурации О, построенное за первые г -1 итераций внешнего цикла, причем о°Р‘ (гк) - «оптимальные» (в оговариваемом далее смысле) плотности тел 0(/), вошедших в представление (2). На шаге г осуществляется перебор (внутренний цикл) всех элементов замощения 0(/'), не вошедших в представление (2). Для каждого из них определяется «оптимальное» значение
плотности о(1) = о°Р‘ (/), доставляющее минимум функционалу
N
ф[о( о]=£
(и(х,) - и(г-1)(х1)У -о(1)и( п(1);х 1)
где
и(г-1)(Х}) = £1о°р (1к)и(П(1к;Х]).
(3)
(4)
к=1
Далее, среди элементов последовательности {г°Р‘ (о}, 1 е {1,2,..., N} \ {г-1, /2,...,^г-1},
определим тот - обозначим его как о °р (/,) -на котором достигнуто минимальное значение функционала (3) (если таких элементов несколько, то берем любой из них). Принимаем О г = О г-1 ^0(/г). Критерием останова ите-
плоскостями рационного
Z — Zл Z — Z -1 Z — Z г
1 2Ь, и в пределах каждой
из полос
I = 1,2,...,Ь -1, осуществлять замощение множества S П1 равновеликими геометри-
ческими телами, одно из которых удовлетворяло бы условию (1).
Собственно алгоритм построения £ -эквивалентного распределения масс представляет собой двухуровневый итерационный
оцесса служит неравенство ф[о°Р‘ (1Г )] <£0 , где £0 - требуемая величина среднеквадратического уклонения наблюденного поля и и поля подобранного «фиктивного» распределения масс по объему О — Ог с кусочно-постоянной плотностью
о(х) = °°р‘(1к), х еО(/'кКк = 1,2,...,п.
Остается добавить, что нулевое приближение О0 = 0 и и(0)(X,).и(X,).
2
Если в рамках заданного замощения не удается «выйти» на распределение масс, обеспечивающих требуемую величину невязки е0, то итерационный процесс повторяется
с другим «подобным» замощением, размеры элементов которого выбираются, как и прежде, в соответствии с условием (1), но при меньшем значении с. При необходимости итерационный процесс может быть выполнен нужное число раз до выхода на требуемое значение невязки. Автоматическая «настройка» алгоритма на нужную структуру замощения позволяет говорить о его адаптивных свойствах.
Несколько замечаний по поводу особенно -стей предлагаемого алгоритма. Фиксированная геометрия локальных «фиктивных» источников модельного поля выводит задачу в разряд линейных, при этом оптимизируется не вектор плотностей
ст = (а(Х),а(2),...,а(п)) (что отвечало бы традиционному подходу к проблеме оптимизации), а последовательно его отдельные ком -поненты ст(/). Разумеется, в последнем случае точность аппроксимации наблюденного поля (при том же числе элементов замощения) несколько теряется. Однако точность, достигнутая при минимизации вектора плотностей , может оказаться выше требуемой (пе-реаппроксимация), а размерность п аппрок-симационной конструкции выше той (т ), что будет установлена по завершении работы предлагаемого алгоритма. Последний является в этом смысле более гибким.
Предложенный алгоритм реализован программно и апробирован на модельном примере. На рис.4 представлено аномальное гравитационное поле, измеренное на резко расчлененном рельефе земной поверхности. Число точек N в числовой модели поля равно 4225. Как показали расчеты, при заданном значении невязки ±0,05 мГал модельное поле обуслов -лено 167 локальными источниками.
Библиографический список
1. Аронов В. И. Обработка на ЭВМ значений аномалий силы тяжести при произвольном рельефе поверхности наблюдений. М.: Недра, 1976. 131 с.
2. Булах Е.Г., Шиншин И.В. Прямые и обратные задачи гравиметрии для совокупности локальных объектов и построение аналитической модели исходного поля // Докл. НАН Украины. 1999. № 1. С. 112-115.
В заключение можно отметить, что новый алгоритм позволяет ускорить процесс аппроксимации геопотенциальных полей для данных большой размерности. Особый интерес представляет использование алгоритма при формировании баз данных в геоинформационных системах [5]. Вместо массива значений наблюденного поля предлагается сохранять аппроксимацию, имеющую существенно меньшую размерность, т.е. перейти от цифровых карт полей к аналитическим аппроксимациям.
Наблюденное поле Рельеф
Рис.4. Иллюстрация работы нового алгоритма построения аппроксимационной конструкции
Аппроксимация обеспечивает возможность восстановления поля в произвольно выбранных точках пространства вне источников, а также выполнение асимптотически оптимальных по точности трансформаций, учитывающих различия в аппликатах ъ точек измерения поля.
3. Бычков С. Г. Методы обработки и интерпретации гравиметрических наблюдений при решении задач нефтегазовой геологии / УрО РАН. Екатеринбург, 2010. 187 с.
4. Долгаль А.С., Пугин А.В. Построение аналитических аппроксимаций геопотенциальных полей с учетом их фрактальной структуры // Докл. РАН. 2006. Т. 410. С. 1152-1155.
5. Долгаль А.С., Симанов А.А. Применение кратномасштабного вейвлет-анализа при аналитических аппроксимациях геопотенциальных по-
лей // Докл. РАН. 2008. Т. 418, № 2. С. 256261.
6. Страхов В.Н. Общая схема и основные итоги развития теории и практики интерпретации потенциальных полей в СССР и России в ХХ веке // Развитие гравиметрии и магнитометрии в ХХ веке: тр. конф.. Москва, 23-26 сентября 1996 г./ ОИФЗ РАН. М., 1997. С. 98-120.
7. Страхов В.Н., Степанова И.Э. Новый информационный базис гравиметрии, магнитометрии и геодезии // Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрических полей: матер. 30-й сессии Международного семинара им. Д.Г. Успенского / ОИФЗ РАН. М., 2003. Ч. I. С. 118123.
Algorithms of Sourcewise Approximation of Geopotential Fields for High Dimension Data
A.F. Sharhimullin, A.S. Dolgal
Perm State University, 614990, Perm, Bukirev st., 15 E-mail: dolgal@mi-perm.ru, art-perm@mail.ru
Some questions concerned with analytic approximations of geopotential fields for a great number of station points are considered. The article includes descriptions of effective algorithms, examples of model and practical interpretative constructions.
Crucially new algorithm of sourcewise structures construction is offered. This algorithm is found on finite element method. Directions of follow-up studies in given district are planned.
Key words: sourcewise approximation, transformation, sources of field, quadtree, wavelet-analysis, finite element method.
Рецензент - доктор технических наук И.А. Санфиров
