Моделирование геологических границ и геопотенциальных полей с использованием быстрого вейвлет-преобразования

УДК 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

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. Сжатие информации о рельефе земной поверхности:

а - исходная ЦММ рельефа

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

Рассмотрим аппроксимации внешних элементов гравитационного и магнитного полей совокупностями эквивалентных источников [1, 2, 7]. Эти аппроксимации называют также «истокообразными аппроксимациями» [14, 15], построением «числовых моделей» [13] или «аналитических моделей» [4] геофизических полей.

Суть такого рода аппроксимационных преобразований заключается в том, что вся информация о наблюденном геофизическом поле U

сохраняется в виде некоторого числа к векторов Р = {р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