CC BY
132
УДК 550.831.016 А.С. Долгаль
Горный институт УрО РАН, Пермь
МОДЕЛИРОВАНИЕ ГЕОЛОГИЧЕСКИХ ГРАНИЦ И ГЕОПОТЕНЦИАЛЬНЫХ ПОЛЕЙ С ИСПОЛЬЗОВАНИЕМ БЫСТРОГО ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ
Вейвлет преобразование - это новое направление в обработке сигналов, бурное развитие которого началось в середине 80-х годов прошлого века. Вейвлеты или «всплески» - функции с компактным носителем, допускающие масштабирование и смещения, образующие полный ортонормированный базис. Иерархическое представление совокупности экспериментальных данных с помощью вейвлетов позволяет описать сигнал произвольной формы в терминах грубого приближения и уточняющих его деталей. Результатом вейвлет-анализа является не только разложение сигнала по частотам (масштабам), но и сведения о временных (пространственных) координатах, на которых эти частоты проявляются [8,9,12].
В геофизике вейвлеты используются, преимущественно, при обработке сейсмических и акустических сигналов. В статье будут рассмотрены возможности применения одного из разновидностей вейвлет-анализа -быстрого вейвлет преобразования (FWT) в качестве эффективного инструмента сжатия информации при кусочно-призматической аппроксимации геологических поверхностей и истокообразной аппроксимации потенциальных геофизических полей.
Быстрое вейвлет-преобразование - FWT (Fast Wavelet Transform), которое иногда также называют кратномасштабным анализом и алгоритмом Малла (Mallat algorithm) [9]. Принцип FWT заключается в том, что для создания «грубого образа» сигнала f(x) служит скейлинг-функция (р(х):
(p(x) = 42Z Икф(2х - к), где к - целые числа; а «уточнение» этого образа
к
происходит с помощью вейвлет-функции \у(х): щ(х) = 42Z gkw(2x - к).
к
Базисные функции (p(x) и цт(х) представляют собой ортонормированные масштабируемые функции с компактным носителем («всплески»), перемещаемые по оси x и легко адаптирующиеся к локальным особенностям сигнала.
Рекурсивное использование процедуры свертки сигнала с коэффициентами hk и gk происходит с уменьшением количества отсчетов в 2 раза при переходе от одного уровня («масштаба») j к другому j+1. Таким образом, происходит отображение сигнала из области его задания в семейство замкнутых вложенных подпространств Vj сVj+1 сVj+2...,
элементами которых являются ортогональные функции (р(х) и щ(х):
f (х) = Z sjn,kPjn ,к + Z Z dj,kVj,k , где sjn,к, dj,к- вейвлет-коэффиЦиенты, jn
к j < jn к
- заданный максимальный уровень разложения сигнала. Аппроксимирующие вейвлет-коэффициенты определяются по итерационной формуле Я]+1 к = Х ^тЯ] 2к+т , Детализирующие коэффициенты - по формуле
т
+1 к = X §тЯ] 2к+т , при этом вся информация о сигнале /(х) сохраняется в
т
} П
наборе коэффициентов ] к и X к (рис. 1).
п’ ] =0 ’
f(x)
d3
Рис. 1. Структура представления сигнала при быстром вейвлет-
преобразовании (заштрихованы коэффициенты, использующиеся при
реконструкции сигнала)
Для приближенного восстановления (синтеза) сигнала
fs(x) ~Zsjn,kPjn,k + Z Zdj/j,к
k j < jn k
вполне допустимо отбросить некоторое количество п0 сравнительно малых коэффициентов dj к, удовлетворяющих условию \djk |<s. При
этом точность восстановления сигнала в метрике L2 будет определяться
выражением s =|| f (x) - fs(x)|| < S^[n0.
Установлено, что вейвлет-коэффициенты существенно отличаются от нуля только вблизи сингулярностей f(x), т.е. вейвлет-ряды обычных функций допускают сильное «разряжение». По этой причине обработка сигнала с помощью FWT позволяет существенно сжать объем информации, отбросить мелкие детали и выделить его наиболее существенные особенности [8].
В качестве базисных функций будем рассматривать функции Хаара (Haar) cp(x) = 9(x)9( 1 — x) и /(x) = 9(x)9(1 — 2x) — 0(2x — 1) 9(1 — x) , /(x) = 9(x)9(1 - 2x) - 9(2x -1)0(1 - x), где 3(x) - функция
Хевисайда (9(x) = 0 при x < 0, 9(x) = 1 при x > 0), а условия на границах имеют вид p(0) = 1,p(1) = 0 и /(0) = 1 , /(0.5) = —1, p(1) = 0. Двумерный базис FWT формируется путем тензорного произведения функций одномерного базиса, в частности для т.н. нестандартного базиса [12, 16] используется единственная скейлинг-функция
pj^ = 2jn pp(2Jn x - к, 2Jn y -1) и три вейвлета
p/h = 2jp/(2jx — k, 2jy — l), /Pki = 2J/p(2J x — k, 2j y — l),
и
цпу]к1 = 2 /(х)^(х)(2] х - к,2 у -1), где к, I - горизонтальный
вертикальный сдвиги, соответственно (рис. 2).
Аппроксимация реального
распределения физических
неоднородностей геологической среды набором тел правильной геометрической формы является необходимым этапом при аналитическом решении прямой задачи гравиразведки [5,6]. Одна из широко используемых аппроксимационных
стратегий базируется на плотном заполнении моделируемых геологических тел прямоугольными призмами, которые соприкасаются боковыми гранями, но не пересекаются друг с другом. В рамках данной стратегии сравнительно легко выполняется геометризация весьма сложных типов природных распределений масс [10,11]. При некотором удалении точки расчета поля от призмы допускается использование приближенных выражений для определения аномального эффекта. Это позволяет существенно ускорить процесс вычислений, практически не ухудшая точность решения прямой задачи [11].
Однако полностью избежать ограничений, связанных с размерностью решаемых задач, при моделировании сложных геологических сред, пока не представляется возможным.
Для решения прямых задач большой размерности крайне важно оптимизировать количество аппроксимирующих прямоугольных призм при сохранении требуемой точности е описания геологических границ т.е.
N ~
V-!~
'=1 I,
Рис. 2. Нестандартные базисные функции Хаара:
А - скейлинг-функция у/ ; вейвлеты: Б - (ру/}к,, В - у/у/1,, Г -у/ф: (светлые участки отвечают
значениям функций +1, темные участки - значениям -1)
добиться выполнения условия
<е, где V- объем моделируемого
геологического тела; Vi - объем ¡-ой призмы, ¡ = 1, 2, ... , N при минимальном количестве призм N [6].
Двумерное с базисными функциями Хаара <рЦ,0, (ру]к ¡, ¥Рк /, ^//
можно рассматривать как кусочно-призматическое представление некоторой геологической границы (например - структурной поверхности, заданной в виде матрицы глубин {Иу} в узлах равномерной сети точек). Количество прямоугольных призм при фиксированной точности е описания геологической поверхности зависит от значений вейвлет-коэффициентов к, большинство которых обычно оказываются пренебрежительно малыми
по модулю. На данной основе можно конструировать весьма технологичные алгоритмы решения прямых задач гравиразведки большой размерности, ускоряющие вычислительный процесс в десятки и сотни раз.
Рассмотрим пример практического использования Б^Г для разряжения сети высот цифровой модели местности (ЦММ). Матрица высот рельефа местности для одной из нефтеперспективных площадей Западного Урала имеет размер 256 строк, 256 столбцов; перепад высот составляет 418 м, при среднем значении 220.6 м и среднеквадратическом отклонении +61.4 м. В результате применения Б'^Г установлено, что с использованием всего лишь 8% исходной информации поверхность рельефа может быть восстановлена с погрешностью +22.7 м. Соответственно, более чем в 10 раз повышается скорость решения прямой задачи гравиразведки от упрощенной ЦММ, сохраняющей все основные особенности рельефа (рис. 3).
а б
Ом 100м 200м 300м 400м 500м 600 м
Рис. 3. Сжатие информации о рельефе земной поверхности:
а - исходная ЦММ рельефа (65536 значений высотных отметок); б - ЦММ, построенная с использованием быстрого вейвлет-преобразования (5171 значение высотных отметок)
Быстрое вейвлет-преобразование для оптимизации кусочнопризматической аппроксимации геологических объектов также может успешно применяться при решении прямых задач в целом ряде других геофизических методов (например - в магниторазведке, в электроразведке методом заряда и т. д.). Особой эффективностью Б'^Г будет обладать при включении его в итерационные алгоритмы метода подбора, базирующиеся на многократном решении прямой задачи от последовательно уточняющихся моделей геологического строения исследуемых территорий.
Рассмотрим аппроксимации внешних элементов гравитационного и магнитного полей совокупностями эквивалентных источников [1, 2, 7]. Эти аппроксимации называют также «истокообразными аппроксимациями» [14, 15], построением «числовых моделей» [13] или «аналитических моделей» [4] геофизических полей.
Суть такого рода аппроксимационных преобразований заключается в том, что вся информация о наблюденном геофизическом поле U(x, у, 2)
сохраняется в виде некоторого числа к векторов Р = {р1, р2,...,рп } параметров источников, создающих модельное поле имод (х, у, 2), практически
эквивалентное полю и(х, у, 2). Определение неизвестных параметров источников выполняется путем решения обратной задачи (ОЗ), которая обычно сводится к минимизации функционала
F(к, Р) = Е [и(X', у', 2') - имод (X', у', 2', к, Р)] на множестве и точек задания '=1
поля. Наибольшие трудности при построении аналитических аппроксимаций поля и(х, у, 2) связаны с большой и сверхбольшой размерностями СЛАУ
[15], т.к. при решении практических задач и ^ 105 -106. Соответственно, для создания эффективных вычислительных алгоритмов, реализующих
«истокообразную аппроксимацию», необходимо использовать минимальное
*
количество к источников поля и (х, у, 2), а также уменьшить число п
параметров, входящих в вектор Р .
Установлено, что пространственное распределение аномалий
гравитационного и магнитного полей обладает приближенной масштабной инвариантностью. Геопотенциальные поля являются мультифракталами, т.к. обладают самоподобной иерархически упорядоченной структурой [3] , которую можно попытаться учитывать в процессе аппроксимации.
Известно, что вейвлет-преобразование, характеризующееся иерархическим базисом, отлично приспособлено для анализа фрактальных и мультифрактальных множеств [8]. Б^Г может эффективно использоваться для декомпозиции исходной задачи аппроксимации поля на несколько последовательно решающихся задач значительно меньшей размерности.
При этом аппроксимации эквивалентными источниками предшествует двумерное Б^Г поля {Дёу}, 1 < ' < т, 1 < у < т, с нестандартными
базисными функциями Хаара. Определяется пороговое значение 8 вейвлет-
коэффициентов йгк I, г = 1, 2, 3, обеспечивающее требуемую точность
восстановления поля е . Затем на самом нижнем уровне ] = ]п проводится
истокообразная аппроксимация поля AgJ т (х, у) = 2 ]п ЕЕ^Ап
к I
(размерность матрицы - т /2п строк, п/2п столбцов) с использованием сеточной модели источников поля: точечные источники (шары) размещаются на глубине I < к < 21 от поверхности Земли (I - расстояние между точками поля). Массы источников М определяются путем решения СЛАУ ОМ = Дg, где О - оператор решения прямой задачи гравиметрии для шара при М = 1. Достаточно хорошая обусловленность СЛАУ обеспечивается указанным выше соотношением ¡/И.
Источники (шары) размещаются только под теми точками задания поля,
г=3
в которых выполняется условие Е ^ГI ^ 38 при данном уровне разложения.
Г=1
На следующих уровнях у < ]п последовательно аппроксимируют эквивалентными источниками «детализирующие» составляющие поля, которые определяются разностью восстановленных по вейвлет -
коэффициентам аномалий и модельного поля Д%^ фрагмента
аппроксимационной конструкции, построенного на предыдущих этапах, причем при размещении источников используется описанный выше принцип. Данный прием позволяет избежать накопления погрешностей
е*,0 < у< Уп, обусловленных собственно истокообразной аппроксимацией на
каждом уровне («масштабе») разложения поля. Циклическое выполнение этих операции в процессе определения параметров эквивалентных источников обеспечивает достижение среднеквадратической погрешности
*
аппроксимации исходного поля эквивалентными источниками е < е.
Рассмотрим пример аппроксимации гравитационного поля Дg шара,
залегающего под хребтообразной формой рельефа, заданного в 64 х 64 = 4096 узлах квадратной сети (т = 64). Диапазон изменения амплитуды поля составляет 0.95 мГал - 17.18 мГал, при среднем значении 5.54 мГал и среднеквадратическом отклонении +3.24 мГал. Максимальный уровень разложения поля при Б'^Г уп = 3. Выбор порогового значения 8 = 0.5
обеспечивает восстановление поля с точностью е=+0.17 мГал с помощью вейвлет-синтеза.
*
Несколько более высокую точность аппроксимации поля е = +0.14 мГал обеспечивает модель, состоящая всего из 342 точечных масс, расположенных на конкордантных рельефу поверхностях с глубинами к=3=237.5 м, к=2 = 118.8 м, к]=1 = 59.3 м. Наличие источников на каждому-ом уровне определяется высокими значениями вейвлет-коэффициентов
г =3
Е ^ I > 1.5, которые в свою очередь, свидетельствуют о наличие
г=1
сингулярностей поля, отчетливо проявленных в соответствующем масштабе (рис. 4).
Рис. 4. Аппроксимация гравитационного поля V с использованием быстрого вейвлет-преобразования
А - синтезированное поле Д^ и аппроксимирующие его источники (58 шаров); Б - синтезированное поле ~ и аппроксимирующие его источники (155 шаров); В -синтезированное поле Д^ 1 и аппроксимирующие его источники (129 шаров); Г - исходное поле Д^ .
1 - эквивалентные источники;
2 - аномалиеобразующий объект
Для многих исследователей, работающих в различных областях науки и инженеров-практиков уже очевидны красота вейвлет-преобразования и его практическая польза. Несомненный интерес представляет использование вейвлетов при решении задач гравиметрии и других геофизических методов. Представленная статья отражает лишь первый, весьма ограниченный, опыт применения этого мощного инструмента применительно к геопотенциальным полям.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Аронов В.И. Обработка на ЭВМ значений аномалий силы тяжести при произвольном рельефе поверхности наблюдений. М.: Недра, 1976. 131 с.
2. Аронов В.И. Методы построения карт геолого-геофизических признаков и геометризации залежей нефти и газа на ЭВМ. М.: Недра, 1990. 300 с.
3. Блох Ю.И. Проблема адекватности интерпретационных моделей в гравиразведке и магниторазведке.//Геофизический вестник. 2004. № 6. с. 10-15.
4. Булах Е.Г., Шиншин И.В. Прямые и обратные задачи гравиметрии для совокупности локальных объектов и построение аналитической модели исходного поля // Докл. НАН Украины. 1999. № 1. С. 112 - 115.
5. Голиздра Г.Я. Основные методы решения прямой задачи на ЭВМ.// Региональная, разведочная и промысловая геофизика. М: ВИЭМС, 1977. 98 с.
6. Гольдшмидт В.И. Оптимизация процесса количественной интерпретации данных гравиразведки. М.: Недра, 1984. 184 с.
7. Долгаль А.С. Аппроксимация геопотенциальных полей эквивалентными источниками при решении практических задач. // Геофизический журнал. 1999. Т. 21. № 4. С. 71-80.
8. Дремин И.М., Иванов О.В., Нечитайло В.А. Вейвлеты и их использование // Успехи физических наук. 2001. Том 171. № 3. С. 465-501.
9. Дьяконов В.П. Вейвлеты. От теории к практике.М: СОЛОН_Р, 2002. 448 с.
10. Матусевич А.В. Объемное моделирование геологических структур на ЭВМ. М.: Недра, 1988. 184 с.
11. Результаты применения моделирования в рудной геофизике в различных районах Сибири. / Под ред. В.С. Моисеева, Г.Г. Ремпеля. М.: Недра, 1988. 219 с.
12. Столниц Э., ДеРоуз Т., Салезин Д. Вейвлеты в компьютерной графике. Пер. с англ. Ижевск: НИЦ «Регулярная и хаотическая динамика», 2002. 272 с.
13. Страхов В.Н. Геофизический «диалект» языка математики. М: ОИФЗ РАН, 2000.
60 с.
14. Страхов В.Н. Основные направления развития теории интерпретации гравиметрических данных в начале XXI века. // Геофизический журнал. 2003. Т. 25. №3. С. 3-8.
15. Страхов В.Н., Степанова И.Э. Новый информационный базис гравиметрии, магнитометрии и геодезии. /Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрических полей. Материалы 30-й сессии Международного семинара им. Д.Г. Успенского. Часть I. М.: ОИФЗ РАН, 2003. С. 118123.
16. Уэлстид. С. Фракталы и вейвлеты для сжатия изображений в действии. М.: Изд-во «Триумф», 2003. 320 с.
© А.С. Долгаль, 2006
