CC BY
193
2014
ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА
Геология
Вып. 4 (25)
ГЕОФИЗИКА, ГЕОФИЗИЧЕСКИЕ МЕТОДЫ ПОИСКОВ ПОЛЕЗНЫХ ИСКОПАЕМЫХ
УДК 550.831
Моделирование геологических объектов и геофизических полей с использованием вейвлетов Хаара
А.С. Долгаль
Горный институт УрО РАН, 614007, Пермь, ул. Сибирская, 78-А E-mail: dolgal@mi-perm.ru
(Статья поступила в редакцию 7 октября 2014 г.)
Рассматривается применение быстрого вейвлет-преобразования с базисными функциями Хаара для моделирования структурных поверхностей и потенциальных геофизических полей, обладающих фрактальными особенностями. Разномасштабное представление экспериментальных данных позволяет существенно снизить затраты при обработке данных большой размерности и повысить качество интерпретационных построений. Представлены алгоритмы кусочно-призматической аппроксимации геологических объектов, предварительной оценки количества эквивалентных источников при аналитической аппроксимации поля, определения намагниченности горных пород верхней части геологического разреза.
Ключевые слова: фрактал, множество, базисная функция, вейвлет Хаара, аппроксимация, поле, аномалия, геологический объект.
Введение
Геологические объекты и создаваемые ими геофизические поля обладают фрактальными особенностями, т.е. характеризуются самоподобной иерархически упорядоченной структурой [1]. В частности, в работе [24] дано определение геологической фрактальности как «дробной размерности и самоподобия геологических объектов в широком диапазоне параметра пространство-время». Следовательно, при анализе геолого-геофизической информации вполне допустимо рассматривать фрактальные множества экспериментально определенных числовых величин, на-
пример, физических характеристик горных пород [12]. Основными свойствами фрактальных множеств являются тонкая структура, допускающая наличие элементов произвольно малых размеров; нерегулярность; самоподобие (включая приближенное или статистическое); дробная (ха-усдорфова) размерность, обычно превышающая топологическую [17]. Эффективный инструмент изучения фрактальных множеств - вейвлет-анализ [10].
Вейвлет-преобразование - это новое направление в обработке сигналов, бурное развитие которого началось в середине 80-х гг. прошлого века. Вейвлеты (wavelets) или «всплески» - функции с ком-
© Долгаль А.С., 2014
пактным носителем, допускающие масштабирование и смещения, образующие полный ортонормированный базис. Иерархическое представление совокупности экспериментальных данных с помощью вейвлетов позволяет описать произвольный сигнал в терминах грубого приближения и уточняющих его деталей. Результатом вейвлет-анализа является не только разложение сигнала по частотам (масштабам), но и сведения о временных (пространственных) координатах, на которых эти частоты проявляются [10].
В геофизике вейвлеты используются преимущественно при обработке сейсмических и акустических сигналов. Имеется также достаточно успешный опыт использования вейвлет-анализа в задачах интерпретации геопотенциальных полей [13, 25, 26, 27]. В статье будут рассмотрены возможности разномасштабного представления данных о геофизических полях и их источниках, реализующиеся с помощью вейвлетов Хаара, позволяющие существенно снизить затраты при обработке данных большой размерности и повысить информативность интерпретационных построений.
Быстрое вейвлет-преобразование
Основная идея вейвлет-анализа заключается в представлении дискретных значений некоторого сигнала в виде ряда
N
f ^) = ^ ciФ {^), где в качестве базисной
i=1
системы функций {Ф i (t) : i = 1,2,..., N}
используются перемещаемые по оси t и масштабируемые однотипные функции с компактным носителем - вейвлеты.
Система {Ф , (t): i = 1,2,..., N } в
большинстве случаев является ортонор-мированной, т.е. отвечает условиям \ 10, при i Ф ]
Ф*, Ф; ={' ■ • 0)
' [1, при i = ]
Вейвлеты Ф ( t ) достаточно хорошо локализованы во временной области и легко адаптируются к локальным особенностям сигналаД). С другой стороны, эти
функции обладают ограниченным спектром частот, поэтому при их масштабировании (растяжении или сжатии по оси t) на основе {Ф i (t): i = 1,2,..., N} можно построить набор из N цифровых фильтров, характеризующихся различными полосами пропускания.
В отличие от спектрального анализа Фурье, где используется амплитудно-частотная область f = J f (t)eiatdt для
представления спектра сигналов, при вейвлет-анализе сигнал преобразуется в область «масштаб-время». Вейвлет-анализ имеет существенные преимущества перед традиционными методами спектрального анализа при обработке нестационарных сигналов [10]. Далее будем рассматривать один из методов вейвлет-анализа - быстрое вейвлет-преобразование (БВП), называемое также алгоритмом Малла (Mal-lat algorithm) [11].
Принцип одномерного БВП заключается в том, что для создания «грубого образа» сигнала f(t) служит скейлинг-функция (p(t):
(t) = 42X h((2t - k) , (2)
k
где k - целые числа, а «уточнение» этого образа происходит с помощью вейвлет-функции у (t) :
W(t) = л/2X gkW(2t - k) . (3)
k
Функции (p(t) и y(t) образуют орто-нормированный базис, рекурсивное использование процедуры свертки сигнала с коэффициентами hk и gk происходит с уменьшением количества отсчетов в 2 раза при переходе от одного уровня («масштаба») j к другому j+1. Таким образом, происходит отображение сигнала из области его задания в семейство замкнутых вложенных подпространств V j с V j +1 с V j + 2 ... , элементами которых являются скейлинг- и вейвлет-функции (p(t) и у (t) :
f (t) = X + XXdj-kWjk , (4)
k j <jn k
где ,к, dj,к - вейвлет-коэффициенты, ]п -
заданный максимальный уровень разложения сигнала. Аппроксимирующие вейвлет-коэффициенты определяются по итерационной формуле
Sj+1,к = X hmSj,2к + т , т
детализирующие коэффициенты - по формуле
d
j+i,k
I г
s
m j ,2 k+m
при этом вся информация о сигнале Д) сохраняется в наборе коэффициентов
,к и Xdhk (рис. 1).
j=0
f(t)
Рис. 1. Структура представления сигнала при быстром вейвлет-преобразовании: синим цветом закрашены коэффициенты, использующиеся при синтезе сигнала
Для приближенного восстановления сигнала /() в формуле (4) целесообразно отбросить некоторое количество п0 сравнительно малых коэффициентов dj к,
удовлетворяющих условию | d] к | < 5 . При
этом точность восстановления сигнала в евклидовой метрике L2 будет определяться выражением ||/^) - /5 ^)|| < . Установлено, что вейвлет-коэффициенты существенно отличаются от нуля только вблизи сингулярностей т.е. вейвлет-ряды (4) обычных, сравнительно гладких функций допускают сильное «разряжение». По этой причине обработка сигнала с помощью БВП позволяет существенно сжать объем информации, отбросить его мелкие детали и выделить наиболее существенные особенности [10].
Базисные функции Хаара
Применительно к интересующему нас классу задач необходимо перейти от временной области задания исходных данных к пространственной: / = / ( х ) или / = / (х, у). Для одномерного случая в
качестве базисных функций будем рассматривать функции Хаара (Haar)
( ( X ) = в (х )в (1 - X )
и
щ(х) = в(х)в(1 - 2х) - в(2х - 1)в(1 - х) , где 3 (х) - функция Хевисайда
(в (х) = 0 при х < 0 , в(х) = 1 при х > 0 ), а условия на границах имеют вид (р(0) = 1, ср(1) = 0 и
у (0) = 1, у (0.5) = - 1, у (1) = 0
[20].
Эти функции представляют собой одно-полярный и разнополярный прямоугольный импульсы, т.е. обладают одновременно свойствами симметрии и ортогональности (в отличие от всех остальных существующих вейвлет-базисов). Первое из этих свойств играет особую роль при параметризации геологических тел и границ. Приведем пример одномерного БВП с функциями Хаара для осциллирующей
функции /(х) = е~°'1х sinx, заданной в 256 точках на отрезке [0, 25.5] оси ОХ с шагом 0.1. Амплитуда/(х) изменяется от -0.63 до 0.85 при математическом ожидании М = 0.04 и среднеквадратическом отклонении S = ±0.31м. БВП при 8 уровнях разложе-
ния сигнала 3 = 8) позволяет получить бенности, при заданной точности и суще-сглаженные представления исходной ственном уменьшении объема информа-функции, сохраняющие ее основные осо- ции (табл. 1, рис. 2).
Таблица 1. Результаты одномерного БВПмодельной функции
Число отбракованных Процент отбракованных Среднеквадратическа,я
погрешно ления f(x)
, j , , j погрешность восстанов-
коэффициентов ujk коэффициентов ujk ^ fAA
183 71.5 0.020
226 88.2 0.050
247 96.4 0.100
т
0 10- 20 30
0 10 20 30
f(x) в
0 10 20 30
Двумерный базис БВП формируется путем тензорного произведения функций одномерного базиса и определяется схемой разложения сигнала, в частности, для т.н. нестандартного базиса [25] используется единственная скейлинг-функция Ро 0 (х, У) = РР(х, У) и три вейвлет-функции
Ф¥к,1 = 23 Р¥(2Зх - кЛ3 у -1\
¥Рк,1 = 2 3 ¥Р(2 3 х - к,2 3 У -
¥¥к1 = 2 :щ(х)щ(х)(23 х - к,23 у - /), где 3
- уровень разложения; к, / - горизонталь-
Рис. 2. Одномерное БВП с базисными функциями Хаара: а - исходная функция f (x) = e ~01 x sin x (256 значений) и результаты ее восстановления после исключения: б - 183, в-226, г-247 вейв-лет-коэффициентов
ный и вертикальный сдвиги соответственно (рис. 3). Нестандартная схема разложения сигнала f = f (x , y ) (последовательное преобразование как строк, так и столбцов матрицы сигнала), определяющая тип базиса, в отличие от стандартной схемы (многократного выполнения одномерного преобразования матрицы значений сигнала сначала по строкам, а затем по столбцам), является более экономичной, т.к. требует в 1.5 раз меньше вычислительных операций [20].
в
Рис. 3. Двумерные функции базиса Хаа-ра для нестандартной схемы разложе-
0 / ния сигнала: а - ; б - ,; в -
¥¥к1> г -уФк1 (желтые участки отвечают значениям функций +1, фиолетовые участки - значениям -1)
Аппроксимация геологических тел и границ
Аппроксимация реального распределения физических неоднородностей геологической среды набором тел правильной геометрической формы является необходимым этапом при аналитическом решении прямой задачи гравиразведки [5]. Одна из широко используемых ап-проксимационных стратегий базируется на плотном заполнении моделируемых геологических тел прямоугольными призмами, которые соприкасаются боковыми гранями, но не пересекаются друг с другом. В рамках данной стратегии сравнительно легко выполняется геометризация весьма сложных типов природных распределений масс [15, 23]. При некотором удалении точки расчета поля от призмы допускается использование приближенных выражений для определения аномального эффекта. Это позволяет существенно ускорить процесс вычислений, практически не ухудшая точность решения прямой задачи [19]. Однако полностью избежать ограничений, связанных с размерностью решаемых задач при моде-
лировании сложных геологических сред, пока не представляется возможным.
При решении прямых задач большой размерности крайне важно оптимизировать количество аппроксимирующих прямоугольных призм при сохранении требуемой точности £ описания геологических границ. т.е. добиться выполнения условия
N ~
V-XV
1=1
<£
(5)
где V - объем моделируемого геологического тела; Vi - объем 1-й призмы, 1 = 1,
2, ... , N при минимальном количестве призм N [4].
Цифровой моделью структурной карты, как правило, является матрица глубин соответствующей геологической границы H = {И/}. Двумерное БВП матрицы H с базисными функциями Хаара
Р0^ (¥к,1, ¥(1,1, ¥¥к,1, благодаря их топологическим особенностям (в т.ч. симметрии), можно рассматривать как кусочно-призматическое представление этой геологической границы. Количество
L
2
Рис. 4. Представление информации о рельефе земной поверхности с помощью БВП: 4 уровня разложения исходной матрицы высотных отметок: для каждого уровня 3 приведены скейлинг
и три вейвлет-коэффициента для базисных функций р|°0, Р¥к
¥Рк1, ¥¥к1 соответственно
прямоугольных призм при фиксированной точности £ описания поверхности (5) зависит от значений вейвлет-коэффициентов dзk, большинство которых обычно оказываются пренебрежительно малыми по модулю. На данной основе можно конструировать весьма технологичные алгоритмы решения прямых
'0,0 > У г к,/■
задач гравиразведки большой размерности, ускоряющие вычислительный процесс в десятки и сотни раз.
Рассмотрим пример практического использования БВП для разряжения сети высот цифровой модели местности (ЦММ). ЦММ для одной из нефтепер-спективных площадей Западного Урала
имеет размер 128 строк, 128 столбцов; перепад высот составляет 393 м при среднем значении 219.9 м и среднеквадратиче-ском отклонении ±61.0 м. На рис. 4 проиллюстрированы результаты ее разложения с помощью БВП для уровней j = 1, 2, 3, 4.
В результате синтеза высот по вейв-лет-коэффициентам при пороговом значении 5 = 40 получена новая ЦММ, отвечающая верхнему уровню j = 5 (рис.5, б). Как очевидно, с использованием лишь около 10% исходной информации поверх-
ность рельефа может быть восстановлена с погрешностью ±13.7 м (т.е. для ее описания теперь достаточно 1696 призм). Можно заметить сохранение всех основных морфологических особенностей рельефа в новой ЦММ (рис.5, б) относительно исходных данных (рис. 5, а). Соответственно приблизительно в 10 раз повышается скорость решения прямой задачи гра-виразведки от упрощенной ЦММ, сохраняющей все основные особенности рельефа (рис.5, б).
100 м
200 м
300 м
-1-
400 м
Рис. 5. Сжатие данных о рельефе земной поверхности с помощью БВП: а - исходная цифровая модель рельефа (16384 значений высотных отметок); б - синтезированная по 1696 вейв-лет-коэффициентам цифровая модель рельефа (исключено ~90% информации)
Рассмотренный выше прием приближенного описания рельефа дневной поверхности может применяться для описания любых структурных поверхностей, а также легко модифицируется применительно к отдельным геологическим телам любой формы. БВП для оптимизации кусочно-призматической аппроксимации геологических объектов также может успешно применяться при решении прямых
задач в целом ряде других геофизических методов (например, в магниторазведке, в электроразведке и т.д.). Особой эффективностью БВП будет обладать при включении его в итерационные алгоритмы метода подбора [2], базирующиеся на многократном решении прямой задачи от последовательно уточняющихся моделей геологического строения исследуемых
разрезов (или 3D-моделей изучаемого объема среды).
Оценка числа элементарных источников при моделировании потенциальных геофизических полей
В настоящее время для моделирования геопотенциальных полей широко используются формальные пространственные распределения эквивалентных источников, создающих близкие наблюденному полю аномальные эффекты (в рамках £-эквивалентности). Определение параметров эквивалентных источников называется также «построением числовых моделей полей» [22] и «истокообразными аппроксимациями» [21]. Для создания эффективных вычислительных алгоритмов, реализующих «истокообразную аппроксимацию», необходимо использовать минимальное количество к элементарных источников поля. Для этого необходимо тем или иным способом согласовать их размещение с морфологическими особенностями наблюденного поля, т.е. провести адаптацию параметров аппроксимацион-ной конструкции применительно к условиям конкретной задачи. Существующие различия в спектральном составе поля в пределах различных подобластей его задания очевидным образом приводят к идее расположения источников на различных уровнях (глубинах). Одной из реализаций этой идеи является аппрокси-мационный алгоритм квадродерева [7], который успешно применяется к данным большой размерности - массивам, содержащим т~ 105-106 значений поля и более. Решение такого рода задач требует значительных затрат машинного времени (до нескольких суток, иногда больше).
Обоснованный выбор параметров для работы алгоритма квадродерева возможен с помощью предварительной процедуры двумерного БВП с функциями Хаара, которая выполняется с целью разложения анализируемого поля на составляющие, отвечающие разноглубинным фрагментам аппроксимационной конструкции. Выде-
ление сингулярностей анализируемого поля позволяет провести приближенную экспресс-оценку качества истокообразной аппроксимации £ при заданном количестве источников к [8].
Поясним работу алгоритма на практическом примере: площадь исследований размером 1600 км расположена в пределах Охотско-Чукотского вулканогенного пояса; колебания высот рельефа земной поверхности составляют 178 - 1533 м; гравитационное поле д, заданное в т=16384 точках с шагом д = 350 м, изменяется в диапазоне от -12.1 до 13.0 мГал (рис. 6, а, б). Рассчитанная аппрок-симационным методом 2-я вертикальная производная поля Vzzz позволила выделить локальные геоплотностные неоднородности, связанные с разрывными нарушениями и не выходящими на дневную поверхность интрузиями диоритов вблизи золоторудного месторождения (рис. 6, в).
Результаты оценок точности аппроксимации с помощью БВП приведены в графах 1-3 табл. 2, фактическая точность - в графах 4,5 табл. 2. Время одного цикла вычислений (для заданного порога 8) составило порядка 10-20 с.
В результате применения описанного выше алгоритма для 3„ =5 построена сеточная модель, состоящая всего из к=3980 источников (шаров), размещенных на глубинах от 350 до 5600 м от земной поверхности. Данная модель обеспечивает точность решения £ = ±0.21 мГал. По мере движения снизу вверх плотность заполнения сеток элементарными источниками заметно уменьшается, при этом отмечается сгущение источников в плане вблизи участков, характеризующихся наиболее сложным характером моделируемой компоненты поля.
Особый интерес представляет собой использование этого алгоритма при формировании баз данных в геоинформационных системах. Вместо массива значений наблюденного поля предлагается сохранять его аппроксимацию (т.е. параметры эквивалентных источников), имеющую существенно меньшую размерность.
а 6 в
Рис. 6. Карты изоаномал гравитационного поля Ag (а); изогипс рельефа местности (б); изолиний 2-й вертикальной производной силы тяжести Vzzz на уровне 2 км (в): 1 - контуры интрузий диоритового состава; 2 - разрывные нарушения; 3 - золоторудное месторождение
Таблица 2. Параметры БВП и истокообразной аппроксимации гравитационного поля золоторудного месторождения
Пороговое значение 5, мГал Процент отброшенных коэффициентов Е ^ , % Погрешность БВП в>, мГал Отношение к/т, % Погрешность истокообразной аппроксимации е, мГал
1 2 3 4 5
0.3 70.4 ±0.11 20.9 ±0.18
0.4 76.7 ±0.14 13.9 ±0.21
0.5 81.0 ±0.16 10.1 ±0.24
0.7 87.0 ±0.22 5.8 ±0.30
1.0 90.9 ±0.28 3.6 ±0.40
Аппроксимация обеспечивает возможность восстановления поля в произвольно выбранных точках пространства вне источников, а также выполнение асимптотически оптимальных по точности трансформаций, учитывающих характер рельефа поверхности наблюдений.
Применение вейвлет-анализа для оценки эффективной намагниченности горных пород верхней части разреза
Результаты измерений напряженности магнитного поля, выполненных в горных областях с интенсивным эффузивным магматизмом, содержат специфические искажения, обусловленные влиянием не-
однородно намагниченного расчлененного рельефа земной поверхности. Эти искажения фиксируются визуально, в виде магнитных аномалий, имеющих ярко выраженную корреляцию с отдельными формами рельефа земной поверхности (хребет, впадина и т.п.). Для учета этих искажений, по аналогии с гравиразведкой, необходим ввод поправок за влияние магнитного рельефа, который в свою очередь требует сведений о пространственном распределении эффективной намагниченности горных пород Jэф = Jэф (х, у) верхней части разреза [16]. Следует отметить, что возможности традиционных методов определения эффективной намагниченности горных пород в естественном залега-
нии, базирующихся на корреляции амплитуд магнитного поля и высот рельефа [14], часто весьма ограниченны.
Предлагается принципиально новый алгоритм оценки эффективной намагниченности верхней части геологического разреза, базирующийся на вейвлет-преобразовании [9]. Алгоритм определения Jэф = Jэф (х, у) основан на сопоставлении между собой относящихся к единому уровню высот «наблюденного» поля и топопоправок 82(вычисленных при
постоянной намагниченности Jo=co„st) с использование БВП с нестандартными двухмерными функциями Хаара (рис. 5). Корреляция (в широком смысле) сингу-лярностей магнитного поля и топопопра-вок предопределяет возможность «локального» определения величины эффективной намагниченности по соотношению амплитуд этих параметров. Очевидно, что в случае пространственного совпадения сингулярностей и 82соотношение
их амплитуд г\ = в каждой точке
будет линейно связано с величиной
Jэф|J0 .
На каждом из к уровней синтеза алгоритма строятся соответствующие приближения пространственного распределения Jkkф . Полученные результаты циклически суммируются друг с другом, таким образом, методом последовательных приближений, строится карта Jэф.
Важную роль в процессе определения Jэф играют пороговые значения, ниже которых соответствующие вейвлет-коэффициенты обнуляются, т.е. принимается отсутствие сингулярностей в данной подобласти. Эти пороговые значения могут быть либо фиксированными («жесткими») во всем пространстве ортонорми-рованных функций, либо плавающими, т.е. определяться исходя из дисперсии этих коэффициентов раздельно для каждого уровня. «Жесткие» пороговые значения для вейвлет-коэффициентов имеют четкий физический смысл и могут быть увязаны с минимально аномальными эф-
фектами в исходном магнитном поле и в поле топопоправок.
Сам процесс БВП не требует значительных вычислительных затрат, поэтому после выполнения нескольких пробных циклов расчета можно выбрать некоторые пороговые значения вейвлет-
коэффициентов (как правило, различные для поля и топопоправок 82), обеспечивающие достаточную точность восстановления сигналов при исключении от 70 до 80-90% коэффициентов их общего числа. Эти «малоамплитудные» коэффициенты отождествляются с помехами и полностью исключаются из процесса синтеза сигнала.
Применение БВП налагает определенные требования на исходные данные -число строк и столбцов матрицы сигнала должно выражаться числом 2Н, где N -целые числа. Чтобы избежать потери информации или избыточной детализации данных аэромагнитной съемки (шаг сети исходных матриц магнитного поля изменять нежелательно), требуется предварительно проводить экстраполяцию данных.
Для магнитных полей со сложной морфологией возможно выполнение предварительной фильтрации - усреднения в скользящих окнах малого размера, нелинейного сглаживания и процедуры отсева «ураганных» значений. Известно, что напряженность магнитного поля и намагниченность горных пород характеризуются логнормальным законом распределения [3]. С целью фильтрации помех выполняется логарифмирование всех результативных значений эффективной намагниченности, которое приводит распределение данных к нормальному, обладающему математическим ожиданием М и дисперсией D. Отбраковка «ураганных» значений может проводится по трем схемам: «сильной», «средней» и «слабой»: соответственно «ураганные» значения параметра заменяются на граничные для интервалов \-yfD ], [- 2л[В] и
Рассмотрим практический пример оценки Jэф по данным аэромагнитной съемки масштаба 1:100 000 для Юнарин-ской площади (около 14 000 кв. км), расположенной в юго-западной части Тунгусской синеклизы. Перепад высот рельефа составляет около 700 м. Сведения о физических свойствах весьма незначительны. Туфогенные образования нижнего триаса практически немагнитны, однако на участках проявления гидротермально-метасоматических изменений их намагниченность возрастает (в среднем течении р. Ниж. Тунгуска, среди практически немагнитных туфогенных образований встречаются породы с намагниченностью до 1.0 А/м). Повышенной намагниченностью обладают туфопесчаники, туфоалев-ролиты - 0.1 - 3.6 А/м, у агломератовых туфов намагниченность достигает 6.0 А/м.
Большую роль в формировании верхней части разреза играют интрузивные траппы нижнего триаса, образующие пластовые тела, реже дайки. Многоярусные пластовые интрузии обуславливают наличие магнитных аномалий, коррелирую-щихся с положительными формами рель-а
ефа. Намагниченность интрузивных траппов от 0.1 до 5.0 А/м. Отмечается наличие обратнонамагниченных траппов, создающих магнитные аномалии отрицательного знака, их намагниченность составляет от 0.1 до 3 А/м.
При выполнении БВП число отбракованных коэффициентов для обоих параметров превысило 80 %, т.е. анализировались ярко выраженные сингулярности полей. На рис. 7 приведены результаты определения Jэф на различных масштабах при фиксированных порогах усечения вейвлет-ряда. Во всех случаях использовались сглаживание и процедура исключения «ураганных» значений. Число точек, в которых отмечалось наличие грубых погрешностей, не превышало 0.02 % общего числа. Таким образом, с помощью БВП удалось получить петромагнитную характеристику верхней части разреза исследуемой площади, что позволило подавить влияние пород трапповой формации на результаты аэромагнитных измерений и выявить ряд аномалий, обусловленных глубинными источниками [2].
б в
Рис. 7. Карты эффективной намагниченности Зэф горных пород верхней части разреза Юна-ринской площади, построенные с использованием БВП: разномасштабные приближения Зф а - уровень 1 (сеть 36.5x36.5 км); б - уровень 3 (сеть 8.25x8.25 км); в - уровень 6 (сеть 1x1 км); красный контур - площадь аэромагнитной съемки масштаба 1:100 000
Следует заметить, что описанный вы- тационных аномалий в редукции Буге ше алгоритм может использоваться и для [18]. определения переменной плотности промежуточного слоя при вычислении грави-
Заключение
Для многих исследователей, работающих в различных областях науки, и инженеров-практиков уже очевидны красота вейвлет-преобразования и его практическая польза. Несомненный интерес представляет использование вейвлетов в области интерпретации данных гравираз-ведки, магниторазведки и других геофизических методов. Представленная статья отражает лишь первый, весьма ограниченный, опыт применения этого мощного инструмента к анализу ограниченных фрактальных множеств, образованных геометрическими параметрами геологических тел и дискретными значениями потенциальных полей.
Неоценимыми преимуществами вейв-лет-преобразования применительно к задачам прикладной геофизики являются: возможность получения независимой разномасштабной (иерархически структурированной) информации; сохранение привязки выявляемых спектральных и иных особенностей сигналов к пространственным координатам; возможность широкого выбора базисных функций, отвечающих свойствам исследуемого сигнала; отсутствие ограничений, налагаемых свойством стационарности (однородности) сигнала, присутствующих в классических методах спектрального анализа; экономичность вычислительных схем. В частности, вейвлеты могут использоваться при решении прогнозно-поисковых геологических задач в сочетании с различными алгоритмами распознавания образов; для оптимизации картографического представления разнородной геофизической информации (т.н. кратномасштабное представление контуров); при томографической псевдоинверсии гравитационного и магнитного полей.
Библиографический список
1. Блох Ю.И. Проблема адекватности интерпретационных моделей в гравиразвед-ке и магниторазведке//Геофизический вестник. 2004. № 6. С. 10-15.
2. Булах Е.Г. Прямые и обратные задачи гравиметрии и магнитометрии. Киев: Наук. думка, 2010. 463 с.
3. Вахромеев Г.С. Основы методологии комплексирования геофизических исследований при поисках рудных месторождений. М.: Недра, 1973. 152 с.
4. Гольдшмидт В.И. Оптимизация процесса количественной интерпретации данных гравиразведки. М.: Недра, 1984. 184 с.
5. Гравиразведка: справочник геофизика /под ред. Е.А. Мудрецовой, К.Е. Весело-ва. 2-е изд., перераб. и доп. М.: Недра, 1990. 607 с.
6. Долгаль А.С. Использование быстрого вейвлет-преобразования при решении прямой задачи гравиразведки // Докл. Академии наук. 2004. Т. 399, № 8. С. 1177-1179.
7. Долгаль А.С., Пугин А.В. Построение аналитических аппроксимаций геопотенциальных полей с учетом их фрактальной структуры // Докл. Академии наук. 2006. Т. 410. С.1152-1155.
8. Долгаль А.С., Симанов А.А. Применение кратномасштабного вейвлет-анализа при аналитических аппроксимациях геопотенциальных полей// Докл. Академии наук. 2008. Т. 418, № 2. С. 256-261.
9. Долгаль А.С., Червоный Н.П. Учет влияния рельефа земной поверхности при аэромагнитных измерениях // Геошформатика. 2008. № 2. С. 58-66.
10. Дремин И.М., Иванов О.В., Нечитайло
B.А. Вейвлеты и их использование // Успехи физических наук. 2001. Т. 171, № 3.
C. 465-501.
11. Дьяконов В.П. Вейвлеты. От теории к практике. М: СОЛОН_Р, 2002. 448 с.
12. Кобрунов А.И. Математические основы теории интерпретации геофизических данных: учеб. пособие. М.: ЦентрЛит-НефтеГаз, 2009. 288 с.
13. Кризский В.Н., Харитонов И.В. Вейвлет-анализ данных геоэлектроразведки постоянным током // Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрических полей: материалы 35-й сессии Меж-дунар. научн. семинара им. Д.Г. Успенского. Сыктывкар, 2007. С. 159-161.
14. Магниторазведка: справочник геофизика / под ред. Б.Е. Никитского, Ю.С. Глебовского. 2-е изд., перераб. и доп. М., 1990. 470 с.
15. Матусевич А.В. Объемное моделирование геологических структур на ЭВМ. М.: Недра, 1988. 184 с.
16. Методические рекомендации по средне-и крупномасштабной магнитной картографии/ НПО «Рудгеофизика». Л., 1990. 84 с.
17. Морозов А.Д. Введение в теорию фракталов / Институт компьютерных исследований. Ижевск, 2002. 160 с.
18. Новикова П.Н. Оценка плотности пород верхней части геологического разреза с использованием вейвлет-преобразования // Геология и полезные ископаемые Западного Урала: матер. регион. науч.-практ. конф / Перм. гос. ун-т. Пермь, 2009. С. 231-235.
19. Ремпель Г.Г. Актуальные вопросы методики введения поправок, связанных с рельефом местности в данные гравираз-ведки и магниторазведки// Изв. АН СССР. Сер. Физика Земли. 1980. № 12. С. 75-89.
20. Старостенко В.И. Устойчивые численные методы в задачах гравиметрии. Киев: Наук. думка, 1978. 227 с.
21. Столниц Э., ДеРоуз Т., Салезин Д. Вейв-леты в компьютерной графике: пер. с англ./ НИЦ «Регулярная и хаотическая динамика». Ижевск, 2002. 272 с.
22. Страхов В.Н. О геофизическом «диалекте» универсального «языка» естественных наук: сб. науч. тр. В 2 т. М.: ОИФЗ РАН, 2002. Т 1. 100 с.
23. Страхов В.Н., Лучицкий А.И. Решение прямой задачи гравиметрии и магнитометрии для некоторых классов распределения масс // Изв. АН СССР. Сер. Физика Земли. 1980. №10. С. 48-64.
24. Сунгатуллин Р.Х. Интегральная геология. Казань, 2006. 142 с.
25. Утемов Э.В., Нургалиев Д.К. "Естественные" вейвлет-преобразования гравиметрических данных: теория и приложения //Физика Земли. 2005. №4. С. 88-96.
26. Христофоров А.В., Абросимова И.С., Христофорова Д.А. Использование вейв-лет-анализа для комплексной интерпретации аномалий теплового и гравитационного полей// Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрических полей: материалы XXXVI сессии Меж-дунар. научн. семинара им. Д.Г. Успенского. Казань, 2009. С. 125.
27. Штокаленко М.Б., Алексеев С Г. Вейв-лет-преобразования с физическим смыслом// Матер. 34-й сессии Междунар. науч. семинара им. Д. Г. Успенского. М., 2007. С. 293-297.
Modeling of Geological Objects and Geophysical Fields Using Haar Wavelets
A.S. Dolgal
Mining Institute of Ural Branch of the Russian Academy of Sciences, 78A Si-birskaya Str., 614007 Perm, Russia. E-mail: dolgal@mi-perm.ru
This article is a presentation of application of the fast wavelet transform with basic Haar functions for modeling the structural surfaces and geophysical fields, characterized by fractal features. The multiscale representation of experimental data allows reducing significantly a cost of the processing of large volume data and improving the interpretation quality. This paper presents the algorithms for sectionally prismatic approximation of geological objects, for preliminary estimation of the number of equivalent sources for the analytical approximation of fields, and for determination of the rock magnetization in the upper part of the geological section. Key words: fractal, aggregate, basic function, Haar wavelet, approximation, field, anomaly, geological object.
References
1. Blokh Yu.I. 2004. Problema adekvatnosti in-terpretatsionnykh modeley v gravirazvedke i magnitorazvedke [Problem of model adequacy in gravimetry and magnitometry]. Geofizicheskiy vestnik, 6: 10-15.
2. Bulakh E.G. 2010. Pryamye I obratnye zada-chi gravimetrii i magnitometrii [Forward and inversion problems of gravimetry and magne-tometry]. Kiev, Naukova Dumka, p. 463.
3. Vakhromeev G.S. 1973. Osnovy metodologii kompleksirovaniya geofizichesrikh issledo-vaniy pri poiskakh rudnykh mestorozhdeniy [Introduction to methodology of integrated geophysical investigation for exploration of mineral deposits]. Moskva, Nedra, p. 152.
4. Goldshmidt V.I. 1984. Optimizatsiya protsessa kolichestvennoy interpretatsii dannykh gravi-razvedki [Optimization of process of quantitative interpretation of the gravity data]. Moskva, Nedra, p. 184.
5. Gravirazvedka: spravochnik geofiizika [Gravimetry: handbook of geophysicist]. Eds. Mudretsova E.A., Veselov K.E. Moskva, Ne-dra, p. 607.
6. Dolgal A.S. 2004. Ispolzovanie bystrogo wavelet-preobrazovaniya pri reshenii prya-moy zadachi gravirazvedki [Usage of fast wavelet transform for solution of gravity forward problem]. Dokl. Akademii nauk, 399 (8): 1177-1179.
7. Dolgal A.S., Pugin A.V. 2006. Postroenie analiticheskikh approksimatsiy geopotentsial-nykh poley s uchetom ikh fraktalnoy struktury [Construction of analytical approximation of geopotential fields taking into account their fractal structure]. Dokl. Akademii nauk, 410: 1152-1155.
8. Dolgal A.S., Simanov A.A. 2008. Primenenie kratnomasshtabnogo wavelet-analiza pri analiticheskikh approksimatsiyakh geopotent-sialnykh poley [Usage of multiple scale wavelet analysis for analytical approximation of geopotential fields]. Dokl. Akademii nauk. 418 (2): 256-261.
9. Dolgal A.S., Chervony N.P. 2008. Uchet vli-yaniya relyefa zemnoy poverkhnosti pri aeromagnitnykh izmereniyakh [Topography corrections for aeromagnetic survey]. Geoin-formatika. 2: 58-66.
10. Dremin I.M., Ivanov O.V., and Nechitaylo V.A. 2001. Wavelety i ikh ispolzovanie [Wavelets and their application]. Uspekhi fizicheskikh nauk. 171 (3): 465-501.
11. Dyakonov V.P. 2002. Wavelety. Ot teorii k praktike [Wavelets. From theory to practice]. Moskva, SOLON_R, p. 448.
12. Kobrunov A.I. 2009. Matematicheskie osnovy teorii interpretatsii geofizicheskikh dannykh [Mathematical basics of theory of geophysical data interpretation]. Moskva, TsentrLit-NefteGaz, p. 288.
13. Krizskiy V.N., Kharitonov I.V. 2007. Wave-let-analiz dannykh geoelektrorazvedki postoyannym tokom [Wavelet analysis of the DC resistivity data]. In Voprosy teorii i prak-tiki geologicheskoy interpretatsii gravitatsion-nykh, magnitnykh I elektricheskikh poley. Materialy 35th sessii Mezhdunar. Nauch. Seminara im. D.G. Uspenskogo, pp. 159-161.
14. Magnitorazvedka: spravochnik geofizika [Magnetometry: handbook of geophysicist]. Eds. Nikitskiy B.E., Glebovskiy Yu.S. Moskva, 1990, p. 470.
15.Matusevich A.V. 1988. Obyomnoe modeliro-vanie geologicheskikh struktur na EVM [3D modeling of the geological structures using computer]. Moskva, Nedra, p. 184.
16. Metodicheskie rekomendatsii po sredne- I krupnomasshtabnoy magnitnoy kartografii [Methodical recommendations on the medium- and large-scale magnetic mapping]. NPO "Rudgeofizika". Leningrad, 1990, p. 84.
17. Morozov A.D. 2002. Vvedenie v teoriyu frak-talov [Introduction to fractals theory]. Izhevsk, Institut kompyuternykh issledovaniy, p. 160.
18. Novikova P.N. 2009. Otsenka plotnosti porod verkhney chasti geologicheskogo razreza s is-polzovaniem wavelet-preobrazovaniya [Assessment of rock density of upper part of geological section using wavelet transform]. In Geologiya i poleznye iskopaemye Zapadnogo Urala. Materialy region. nauch,-prakt. konf. Perm, Perm State Univ., 2009, pp. 231-235.
19. Rempel G.G. 1980. Aktualnye voprosy meto-diki vvedeniya popravok, svyazannykh s re-lyefom mestnosti v dannye gravirazvedki i magnitorazvedki [Actual problems of methodology of gravimetry and magnetometry corrections related to the surface relief]. Izv. AN USSR. Seria Fizika Zemli, 12: 75-89.
20. Starostenko V.I. 1978. Ustoychivye chislen-nye metody v zadachakh gravimetrii [Stable numerical methods in the gravimetry problems]. Kiev, Naukova Dumka, p. 227.
21. Stolnitz E.J., DeRous T.D., and Salesin D.H. 1996. Wavelets for Computer Graphics:
Theory and Applications. Morgan Kaufmann Publishers, San Francisco.
22. Strakhov V.N. 2002. O geofizicheskom "dialekte" universalnogo "yazyka" estestven-nykh nauk [About the geophysical "dialect" of universal language of natural sciences]. T. 1. Moskva, OIFZ RAN, p. 100.
23. Strakhov V.N., Luchitskiy A.I. 1980. Reshenie pryamoy zadachi gravimetrii i magnitometrii dlya nekotorykh klassov raspredeleniya mass [Solution of forward problem of gravimetry and magnetometry for some examples of mass destribution]. Izv. AN USSR, Seria Fizika Zemli, 10: 48-64.
24. Sungatullin R.Kh. 2006. Integralnaya ge-ologiya [Integral geology]. Kazan, p. 142.
25. Utemov E.V., Nurgaliev D.K. 2005. "Estest-vennye" wavelet-preobrazovaniya gravimet-richeskikh dannykh: teoriya i prilozheniya ["Natural" wavelet transforms of gravity data:
theory and application]. Fizika Zemli, 4: 88-96.
26. Khristoforov A.V., Abrosimova I.S., and Khristoforova D.A. 2009. Ispolzovanie wave-let-analiza dlya kompleksnoy interpretatsii anomaliy teplovogo i gravitatsionnogo poley [Application of wavelet analysis for integrated interpretation of anomalies of thermal and gravitational fields]. In Voprosy teorii i prak-tiki geologicheskoy interpretatsii gravitatsion-nykh, magnitnykh I elektricheskikh poley. Materialy XXXVI Sessii Mezhdunar. nauchn. seminara im. D.G. Uspenskogo. Kazan, p. 125.
27. Shtokalenko M.B., Alekseev S.G. 2007. Wavelet-preobrazovaniya s fizicheskim smys-lom [Wavelet transforms with Physical meaning]. Materialy 34th sessii Mezhdunar. nauchn. seminara im. D.G. Uspenskogo. Moskva, pp. 293-297.
