Трансформации потенциальных полей на основе непрерывного вейвлет-преобразования Текст научной статьи по специальности «Физика»

УДК 550.83.016

К.М. Кузнецов1, И.В. Оболенский2, А.А. Булычев3

ТРАНСФОРМАЦИИ ПОТЕНЦИАЛЬНЫХ ПОЛЕЙ НА ОСНОВЕ НЕПРЕРЫВНОГО ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ4

Рассмотрена возможность применения непрерывного вейвлет-анализа на основе ядра Пуассона для обработки и интерпретации данных грави- и магниторазведки. По вейвлет-спектру W(x, h) возможна реконструкция исходного сигнала g(x). Также по вейвлет-спектру возможно построение эквивалентного распределения масс и намагничен-ностей, создающих исходное поле g(x). Указанные возможности позволяют использовать вейвлет-преобразования для решения таких задач, как фильтрация исходных сигналов, продолжение поля и вычисление высших производных полей в верхнем и нижнем полупространствах и редукция магнитных аномалий к полюсу.

Ключевые слова: непрерывные вейвлет-преобразования, потенциальные поля, продолжение поля, вычисление высших производных.

Some possibility of continuous wavelet transformation based of Poisson kernel for applying for processing and interpretation of gravity and magnetic data are considered in this article. Reconstruction of original signal g(x) by wavetet spectrum W(x, h) is possible. In addition, it is possible to build equivalent distribution of masses and intensity of magnetization, which create original field g(x). Mentioned opportunities allow using wavelet transformations for solving such problems like filtration original signals, upward and downward continuation, field derivative calculation and reduction to pole of magnetic anomalies.

Keywords: continuous wavelet transformation, potential fields, field continuation, derivative calculation.

Введение. При обработке и интерпретации данных гравиразведки и магниторазведки применяются различные трансформации исходных полей. В последние 15—20 лет активно развивается новое направление в задачах анализа сигналов, в том числе геофизических, основанное на вейвлет-преобразованиях [Астафьева, 1996; Никитин, Петров, 2010; Могеаи й а1., 1997]. Однако возможности вейвлет-анализа применительно к аномальным гравитационным и магнитным полям еще недостаточно изучены. Цель работы — развитие этого метода применительно к решению «классических» задач преобразований потенциальных полей, таких, как пересчет поля в верхнее и нижнее полупространство, редукция аномального магнитного поля к полюсу, фильтрация и ряд других, описываемых интегралом типа свертки.

Основы теории вейвлет-преобразования. Одна из разновидностей вейвлет-анализа основана на непрерывном вейвлет-преобразовании. Оно осуществляется путем свертки анализируемой одно-параметрической функции с двухпараметриче-ской функцией уИ,х(%) [Астафьева, 1996; Добеши, 2001; Короновский, 2003; Никитин, 2010; Юдин, 2001; Могеаи е! а1., 1997]:

1¥{х,К)= \ gQ:¡)xfx,h{k)d^ (1)

—оо

Черта над функцией ух,И(%) означает ее сопряжение. Функцию Щх, И) часто называют вейвлет-спектром функции g(%). Сама функция ¥х,И(^) получается из базисного (материнского) вейвлета у0(%):

Параметр И — масштаб вейвлетного преобразования, его значение меняется от нуля до бесконечности (ИеЯ+); х (хеЯ) — параметр сдвига, определяющий положение вейвлета на оси О% (оси Ох). Коэффициент

( ) вводится с целью нормировки этой функции.

Функция у0(%) должна удовлетворять определенным требованиям, в частности:

оо оо

= (3)

—оо —оо

Первое условие означает, что среднее значение функции у0(%) и порождаемых ею вейвлетов ух,И(%)

1 Московский государственный университет имени М.В. Ломоносова, геологический факультет, кафедра геофизических методов исследования земной коры, аспирант; e-mail: kirillkuz90@gmail.com

2 Московский государственный университет имени М.В. Ломоносова, геологический факультет, кафедра геофизических методов исследования земной коры, аспирант; e-mail: ivan-obolensky@yandex.ru

3 Московский государственный университет имени М.В. Ломоносова, геологический факультет, кафедра геофизических методов исследования земной коры, заведующий кафедрой, докт. физ.-мат. н.; e-mail: aabul@geophys.geol.msu.ru

4 Работа выполнена при финансовой поддержке РФФИ (грант № 13-05-00143).

должно быть равно нулю, второе — функция ^о(^) ограничена, а соответственно и функции ухИ(^). Кроме того, значения функции у0(^) должны быстро стремиться к нулевым при стремлении аргумента к бесконечности (^±о>). Вид вейвлета зависит от поставленной задачи и может быть несимметричным.

В случае если исходная функция имеет среднее значение, равное нулю, и ее разложение выполнялось с помощью вейвлетов, удовлетворяющих условиям (2, 3), то возможно осуществить обратное непрерывное вейвлет-преобразование, т.е. по функции Щх, И) восстановить (реконструировать) исходный сигнал g(^):

\

, (4)

где — константа, определяемая функцией у0:

С¥ = |

|¥о(«>)|' со

ч/(0,

где фо(со) — Фурье-образ функции у0(х):

оо

у0(а>) = | ~Уо(х)е~1<ах(1х.

(5)

(6)

.

^ —оо '

(7)

1 ь2 И2-£-х)2

Я [(^-х)2 + А2]2

ч*о(£) =

2$

2(Л-х)к

(10)

Вейвлеты и Ух,к (£) соответствуют вычислению 1-й вертикальной производной, а вейвлеты ¥о (£) и — 1-й горизонтальной производной

исходной функции. Очевидно, что при заданном параметре И функция Ж(х, И), рассчитанная с помощью вейвлетов (9, 10), связана с первой частной производной поля g(^) на высоте соотношением

ь/А,

я к2

.

(11)

Таким образом, на основании соотношений (1, 4) могут быть осуществлены различные линейные трансформации исходных полей. Заметим, что вычисление вейвлет-спектра Щх, И) представляет собой при заданном параметре И линейную трансформацию исходного поля g(^) с ядром ухИ(^), описываемую интегралом типа свертки:

Здесь под функцией g1 понимаются функции gx или gz, соответствующие первым частным производным функции g по параметру х или 2.

Аналогично (9, 10) могут быть представлены вейвлеты Пуассона более высоких порядков.

При задаче реконструкции сигнала необходимо знание значений коэффициента Су, которые можно определить на основании соотношения (5). С учетом того, что рассматриваемые базисные вейвлеты разных порядков определяются через частные производные ядра Пуассона (9, 10), то их Фурье-образ имеет вид

ф^'"12 (ш) = п(т)" |ю| е~

(12)

Соответственно значения коэффициента Су будут представлены следующим интегралом:

2

0

ю

ч/со =

Вейвлеты на основе ядра Пуассона. При анализе аномальных гравитационных и магнитных полей логично использовать вейвлеты, построенные на основе ядра Пуассона, как его частные производные:

. (8)

д х дтг\х1 + г1)

Полученные таким образом вейвлеты будут удовлетворять необходимым условиям (3). Так, вейвлеты, соответствующие вычислению первых частных производных поля (вейвлеты 1-го порядка) на высоте И, будут иметь вид [Оболенский, 2011б]:

1Ч2

= 1

ж2(й2(т+п)е~2а

(¡(0 =

со

со

2ЛГ-1„-2 со

(1(0,

(Ы=т+н).

(13)

Аналитическое решение этого интеграла имеет

вид

-

1 2К-\2Ы-2

1 2{2Ы-\)\

. (14)

2 22м

2 2 2

Так, например, для вейвлетов 1-го порядка Сщ=л2/4, для вейвлетов 2-го порядка

_ 21Г32П 23

^ 2^2 2 2

для вейвлетов 3-го порядка

с =„21(543 2 П 215 ¥ 2^22222/ 8

и т.д.

Комплексный вейвлет Пуассона. Для анализа двумерных потенциальных полей в работе [Оболенский, 2011,] был предложен вейвлет, построенный на комплексном представлении ядра Пуассона:

, (15)

С-5

где а = % + ¡С, — комплексная координата положения двумерного точечного источника, который создает эффект, а я = х + гг — комплексная координата точки наблюдения (г — мнимая единица, ось Ох направлена направо, ось Ог — вверх). Производные этой функции по параметру я имеют вид

где т — порядок производной. Основываясь на последнем соотношении, можно записать вейвлет [Оболенский, 2011б]:

1Г<"> =

г.т+1

■т+1

т\

7И+1 '

1т+1т\

ип

>) =ч__

(17)

„—|со|А

С(со) =

-2пе

. о, С(от)(ю) = 0'(о)тС(ю).

ю>0 ю<0'

(18)

О

кт+1

л/й

ю<0'

ш™2тее-Н*; 0

ю>0 со<0

(19)

торое эквивалентно глубине. Предположим, что значение Н1 связано с шагом задания исходного поля (Иг = ¡Ах). Для быстрого вычисления свертки воспользуемся свойством, заключающимся в том, что спектр дискретной функции, получаемой в результате вычисления свертки, равен произведению дискретных спектров сворачиваемых функций. Тогда алгоритм расчета вейвлет-спектра будет включать следующие шаги:

действительная часть которого является четной, а мнимая — нечетной функцией.

Для функции С(х) согласно выражению (6) можно определить спектр Фурье. Поскольку функция С(х) аналитическая, то спектр самой функции и спектры производных имеют вид

В этом случае спектральные характеристики вейвлетов (17) можно представить в виде

^(»0 _ [ют2л;е~НА; ю>0

Отметим, что действительные части этих функций отличны от нуля лишь при положительных ю, а их мнимые части равны нулю.

Алгоритм численной реализации вейвлет-пре-образования. На рис. 1 представлена общая блок-схема расчета вейвлет-спектра исходного сигнала. Будем полагать, что анализируемая функция представлена N дискретными значениями с постоянным шагом Ах между ними. Согласно формуле (1) вейвлет-спектр — результат свертки входного сигнала и вейвлета на каждом значении И, ко-

Рис. 1. Блок-схема алгоритма вычисления вейвлет-спектра

1) расчет дискретного спектра исходного сигнала на основе алгоритма быстрого преобразования Фурье (БПФ по основанию 2). Применение алгоритма БПФ предполагает, что ^2т (т=1, 2, 3, ...);

2) вычисление для каждого значения И соответствующей частотной характеристики вейвлет-преобразования на основе соотношений

(19), юк = (2n/NАx)k, . С учетом свойств

дискретного преобразования Фурье частотная характеристика вейвлет-преобразования рассчитывается для значений частот к = 0,Ж/2, значения при ^N/2 полагаются равными нулю;

3) перемножение полученной частотной характеристики с сопряженными значениями дискретного спектра анализируемого сигнала;

4) вычисление обратного дискретного преобразования Фурье от полученного дискретного спектра, в результате чего образуется комплексная функция при заданном параметре к, действительная часть которой соответствует свертке исходной функции с четным вейвлетом, мнимая — с нечетным.

Трансформации потенциальных полей. На

основании уравнения (4) по полученному вейвлет-спектру Щ(х, к) возможна реконструкция исходного сигнала g(Q. Это в свою очередь позволяет использовать вейвлет-преобразования для решения таких «классических» задач, как фильтрация исходных сигналов, редуцирование двухмерных магнитных аномалий к полюсу, продолжение поля и вычисление высших производных полей в верхнем и нижнем полупространствах.

Остановимся на этих вопросах подробнее.

Построение эквивалентного распределения масс и намагниченности, создающих исходное поле g(^). Пусть функция g(x) представляет собой аномальное гравитационное поле, создаваемое распределением масс и(х, к) в нижней полуплоскости:

§(х) = 2 сЦъ&к)

А

- х)2 + к

(20)

где G — гравитационная постоянная.

Определим функцию ст(х, к), связанную с полем g(x) соотношением:

1 Ъ2 - г>2

, (21)

** (&-х)2+к2)

д(ю,/г) = #(ю)(|со|е-шА).

(22)

Отметим, что полученное плотностное распределение не отображает истинного распределения масс и(х, к), а эквивалентно ему и не несет явного геологического смысла.

Функция ст(х, к) связана с функцией Щ(х, к), рассчитанной через вейвлет (9), соотношением

(11)

ж к2

,

(24)

откуда следует, что эквивалентное распределение плотности 5(х, к), создающее поле g(x), связано с вейвлет-преобразованием следующим образом:

Та

5(х,А) = -^—<з(х,И) = —^

71(7 К Сг А

Щх,К). (25)

Для примера рассмотрена модель, состоящая из двух точечных источников на глубине 100 условных единиц (у.е.), равных расчетному шагу (1х вдоль профиля, с одинаковой по модулю эффективной массой и разной по знакам. На рис. 2 представлен результат вейвлет-преобразования функции g(x), полученной как гравитационный эффект от модели вейвлетом (9), и соответствующее эквивалентное распределение плотности 8(х, к), рассчитанной по формуле (25).

Предположим, что функция g(x) соответствует полю ^(х), создаваемому плотностным распределением масс и(х, к):

А2 - - х)2

§(х) = 2в\ к)

'о (&-х)2+к2)

й^Й?А. (26)

т.е. эта функция совпадает с первой вертикальной частной производной поля g на высоте к. Используя представление (6) для функций g(x), ст(х, к) и ядра преобразования, это выражение в частотной области будет иметь вид

В то же время согласно (4) поле g(í) можно реконструировать по результатам вейвлет-преобра-зования Щ(х, к), полученного вейвлетом (9).

Осуществим следующие преобразования: / „ л

йк =

Придадим функции ст(х, к) смысл плотности и определим поле/(х), создаваемое ею. В частотной области для функции /(ю) можно записать:

оо

/(со) = 2с|д(со,А)(71е-0)Л)^А = о

оо

= 2Сг | ¿(оо)(| со| и А )(7ге~ ш л }йк = о

оо

= 2кв\(а\ #(со)|е-2в) нйк = яС^(ю). (23) о

Из полученного результата следует, что функция 5(х,к) = (1Дб)ст(х,к) описывает такое распределение масс, которое создает исходное поле g(x) [Кобрунов, Варфоломеев, 1981].

т О V

2 Ь2-£-х)

= 2вх

оо оо

*Я

0^2вСук24к

= 2(?| | 8(х,к)

ЩЬЛ)

2 >

(&-х)2 + к2)2

ь2 А2 - - х)2

((%-х)2+к2)2

(27)

где

6(х,А) =

2вСу4к

Как и в предыдущем случае, полученное плот-ностное распределение эквивалентно истинному и(х, к). На рис. 3 показано исходное поле от той же модели, которая использована на рис. 2, и эквивалентное распределение плотности, построенное по формуле (27) в предположении, что функция g(x) соответствует полю У2Г

Рис 2. Результат прямого и обратного вейвлет-преобразования четным вейвлетом 1-го порядка (А); пересчет вейвлет-спектра в

значения плотности и рассчитанное от них поле Vz (Б)

Пусть функция g(x) соответствует вертикальной составляющей индукции аномального магнитного поля Z(x). Воздействуем на эту функцию четным вейвлетом Пуассона 1-го порядка (9). На основании соотношения (4) полученной функции Ж(х, И) можно поставить в соответствие функцию распределения вертикальной намагниченности 1(х,И) в нижнем полупространстве, создающей поле Z(x)

Редукция магнитных аномалий к полюсу. Пусть аномальное магнитное поле создается двухмерным изолированным источником, расположенным ниже оси Ох, с поперечным сечением D и с постоянной намагниченностью. Определим функцию У2(х, г) следующим образом:

" С"*

Vz{x,z) = l\

-d^dC,. (30)

I{x,h) =

1

2 C^yfh

W{x,h\

(28)

поскольку (в системе СИ) g(x) = Z(x) =

= ^2 f f/CE h) *2-(S-*>2

4* Va Ж-*?

dhdÇ,. (29)

где ц0 — магнитная проницаемость вакуума (ц0 = = 4п-10-7 Гн/м).

Полученные результаты можно распространить на вейвлет .

Тогда на оси Ох для вертикальной компоненты магнитной индукции Z(x), создаваемой этим источником, справедливо следующее соотношение (с точностью до постоянного сомножителя (ц0/4п) в системе СИ):

ад(*)+%(*), (31)

где 1х и — компоненты вектора намагниченности вдоль осей Ох и Ог соответственно, а Ух2 и Ук — частные производные от функции У2 по параметрам х и г. Введем следующее понятие эффективной намагниченности объекта I: 1х = /ео8ф,

Рис. 3. Результаты прямого и обратного вейвлет-преобразования четным вейвлетом 1-го порядка (А); пересчет вейвлет-спектра

в плотности и рассчитанное от них поле (Б)

1г = 1зшф, ф — угол ее наклонения. В частотной области выражение (31) можно представить в виде

Z(to) = Ixi(ùVz (со) + Iz(ûVz (to) = = /(sin ср + г cos ф)со Vz (со),

(32)

где Fz((ö) — Фурье-образ функции Vz(x).

Воздействуем на функцию Z(x) вейвлетом 1-го порядка следующего вида:

Va,* Й) = sin <p yzx>h (Ç) + cos ф \fxxJt (£). (33)

(1, 7), для вейвлет-спектра W(x,h) в частотной области можем записать следующее:

- « h2 W{ю, h) = I (sin ф + i cos q>)(ùV2 (to) • к —¡= x

VA

h2

(35)

л/h

Полученной функции Ж(х, И) можно придать смысл интенсивности намагниченности, создающей исходное поле Z(x). Если теперь к С учетом того, что функции и ФункДип Ж(х, И) применить обратное вейвлет-

определяются соотношениями (9, 10), в частотной области вейвлет (33) будет представлен в виде

Va,* (ю) = мп Ф

n—j= toe

4h

-wh

+ СОвф

h2 . n —¡=i(oe

4h

-юА

= %—j=(ùe (sin ф + г совф).

у! h

(34)

Осуществив непрерывное вейвлет-преобра-

преобразование (4) с вейвлетом то по-

лученный результат будет соответствовать вертикальной составляющей аномального магнитного поля ZB(x) в предположении, что намагниченность имеет только вертикальную компоненту. Если же воспользоваться вейвлетом V*,a(£), то результат будет соответствовать горизонтальной компоненте аномального поля XB(x).

Для примера на рис. 4 показан рассчитанный магнитный эффект от точечного источника на

зование функции Z(x) с вейвлетом (33) согласно глубине 100 м при наклонении 60°

Фильтрация полей. Поскольку по значениям вейвлет-спектра \¥(х, И) на основе соотношения (4) может быть восстановлена исходная функция g(^), фильтрация поля (сигнала) в этом случае будет основана на том, что некоторые из значений функции Ж(х,И) обнуляются или изменяются их амплитудные значения. При решении этой задачи можно использовать не только вейвлеты Пуассона, но и другие типы вейвлетов [Оболенский, 20116], которые в ряде случаев могут оказаться предпочтительнее.

В качестве примера рассмотрим результаты вейвлет-фильтрации модельного сигнала. В качестве исходного сигнала рассмотрим поле, создаваемое той же моделью, что и на рис. 1, осложненное случайной помехой. На рис. 5 представлены результаты восстановления сигнала для случая, когда значения функции Ж(х, И) были обнулены для И = 1...20. Фильтрация проводилась вейвлетом Пуассона 3-го порядка (10).

Продолжение и вычисление высших производных поля в верхнем полупространстве. Решение этой задачи тесно связано с задачей определения эффективной плотности или намагниченности по результатам непрерывного вейвлет-преобразования исходного поля. Определив распределение эффективных источников в нижнем полупространстве, согласно (25, 27, 28) можно рассчитать значения поля или его производные на заданной высоте. Формально можно осуществить пересчет полученных значений функции Ж(х, И) в значения функции Ж(х, И), у которой значения будут равны нулю при И<2, где г — высота пересчета поля, а при И>г будут соответствовать коэффициентам вейвлет-преобразования. Однако в данном случае эта процедура излишняя.

На рис. 6 представлены результаты пересчета поля на высоту Н=100 у.е. длины. На этом рисунке показано рассчитанное точное поле на данном уровне и поле, полученное на основе вейвлет-трансформации в верхнее полупространство.

Рис. 5. Результат фильтрации сигнала четным вейвлетом Пуассона 3-го

порядка

Рис. 4. Результаты пересчета к полюсу (А) и вейвлет-спектр от исходного поля (Б)

Продолжение и вычисление высших производных поля в нижнем полупространстве. Решение этой задачи, как и в предыдущем случае, основано на том, что по результатам вейвлет-преобразования с помощью вейвлетов Пуассона Ж(х, И) можно определить эквивалентное распределение плотности 8(х, И) в нижнем полупространстве, создающее исходное поле. Затем, обнулив значения функции 8(х, И) до необходимого уровня пересчета, можно рассчитать прямой гравитационный (магнитный) эффект на заданном уровне. При решении этой задачи исходный сигнал предпочтительнее рассматривать как поле, соответствующее более высоким частным производным гравитационного потенциала, чем поле V, т.е., например, как поле Vzzz. В этом случае эффективные массы окажутся расположенными на большей глубине по сравнению со случаем применения вейвлетов 1-го порядка. На рис. 7 представлены результаты такого пересчета на глубину И=50 у.е. глубины.

Заключение. Таким образом, показано, что при анализе двухмерных потенциаль-

-3000- ^

-положительные изолинии

-3500— - нулевая изолиния

---отрицательные изолинии

-40000 500 1000 1500 2000 2500 3000 3500 4000

х, У- е.

Рис. 6. Результаты пересчета вверх на 100 у.е. глубины с использованием четного вейвлета 1-го порядка

ных полей можно успешно использовать алгоритмы, построенные на непрерывном вейвлет-преобразовании. Пересчеты в верхнее и нижнее полупространства основаны на расчете эквивалентного распределения значений плотности или намагниченности и решении прямой задачи от них на заданном уровне. Редуцирование магнит-

СПИСОК ЛИТЕРАТУРЫ

Астафьева Н.М. Вейвлет-анализ: основы теории и примеры применения // УФН. 1996. Т. 166, № 11. С. 1145-1170.

Добеши И. Десять лекций по вейвлетам. Ижевск: НИЦ «Регулярная и хаотическая динамика», 2001. 464 с.

Дьяконов В.П. Вейвлеты. От теории к практике. М.: СОЛОН-Р, 2002. 448 с.

Кобрунов А.И., Варфоломеев В.А. Об одном методе е-эквивалентных перераспределений и его использовании при интерпретации гравитационных полей // Изв. АН СССР. Физика Земли. 1981. № 10. С. 25-44.

Короновский А.А., Храмов А.Е. Непрерывный вейвлет-анализ и его приложения. М.: Физматлит, 2003. 176 с.

Никитин А.А., Петров А.В. Теоретические основы обработки геофизической информации: Учеб. пособие. М.: ООО «ЦИТвП», 2010. 114 с.

Оболенский И.В., Булычев А.А. Непрерывное вейвлет-преобразование гравиметрических и маг-

X, у. е.

-положительные изолинии

— - нулевая изолиния --- отрицательные изолинии

Рис. 7. Результаты пересчета вниз на 50 у.е. глубины с использованием четного вейвлета 3-го порядк

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

нитометрических данных // Геофизика. 2011а. № 3. С. 48-56.

Оболенский И.В., Булычев А.А. Применение комплексного непрерывного вейвлет-преобразования Пуассона для определения источников аномалий потенциальных полей // Геофизические исследования. 2011б. Т. 12, № 3. С. 5-21.

Утёмов Э.В., Нургалиев Д.К. «Естественные» вейвлет-преобразования гравиметрических данных: теория и приложения // Изв. РАН. Физика Земли. 2005. № 4. С. 88-96.

Юдин М.Н., Фарков Ю.А., Филатов Д.М. Введение в вейвлет-анализ: Уч.-практ. пособие для системы дистанционного образования. М.: МГРА, 2001. 72 с.

Moreau F., Gibert D., Holschneider M., Saracco G. Wavelet analysis of potential fields // Inverse Problems. 1997. Vol. 13. P. 165-178.

Поступила в редакцию 12.05.2015