CC BY
110
УДК 550.83.016
ПРИМЕНЕНИЕ ЭМПИРИЧЕСКОЙ МОДОВОЙ ДЕКОМПОЗИЦИИ ПРИ ОБРАБОТКЕ ГЕОФИЗИЧЕСКИХ ДАННЫХ
Долгаль Александр Сергеевич1,
Христенко Людмила Анатольевна1,
1 Учреждение Российской академии наук «Горный институт Уральского отделения РАН», Россия, 620007, г. Пермь, ул. Сибирская, 78-А.
Актуальность исследований обусловлена целесообразностью использования передовых математических методов при обработке результатов полевых геофизических наблюдений.
Цель работы: повышение эффективности методов прикладной геофизики за счет адаптивного выделения информативных составляющих физических полей с использованием эмпирической модовой декомпозиции (Empirical Mode Decomposition - EMD). Метод исследований: декомпозиция профильных геофизических данных, при которой в ходе итерационного вычислительного процесса адаптивно выделяются ортогональные разночастотные компоненты сигнала (EMD), которые называются эмпирическими модовыми функциями (IMF). Метод предназначен для частотного представления нестационарных сигналов, его отличительной особенностью является отсутствие априорно заданного функционального базиса разложения. Предлагается оригинальный алгоритм, в котором для построения функций, огибающих экстремумы анализируемого сигнала, используется истокообразная аппроксимация, а остаток при разложении может отождествляться с фоновой составляющей поля. Результаты. Алгоритм EMD программно реализован и апробирован на практических материалах при обработке данных грави-разведки, магниторазведки и электроразведки. Установлена возможность применения метода с целью выделения информативной компоненты в результатах электропрофилирования (Пермский край); независимой оценки качества высокоточной гравиметрической съемки (Западный Саян); определения составляющей магнитного поля, обусловленной влиянием резко расчлененного рельефа земной поверхности в условиях развития эффузивных траппов (Норильский район). В последнем случае используется построение серии регрессионных зависимостей между разночастотными составляющими магнитного поля и высот. Сделан вывод о целесообразности использования метода эмпирической модовой декомпозиции в области прикладной геофизики. Перспективы применения EMD авторы видят в возможности анализа материалов геофизического мониторинга разрабатывающихся месторождений полезных ископаемых.
Ключевые слова:
Сигнал, эмпирическая модовая декомпозиция, алгоритм, аппроксимация, составляющая поля, гравиразведка, магниторазведка, электроразведка.
Введение
Метод эмпирической модовой декомпозиции (Empirical Mode Decomposition - EMD) был предложен Норденом Хуангом в 1995 г. и первоначально использовался при изучение поверхностных волн тайфунов. В 1998 г. метод был обобщен применительно к анализу произвольных временных рядов [1]. Метод является важнейшей составляющей преобразования Гильберта-Хуанга (Huang-Hilbert Transform - HHT), получившего в дальнейшем широкое применение в различных областях науки и техники [2-5 и др.], наряду с преобразованием Фурье [6] и вейвлет-анализом [7]. В области геофизики известны примеры успешного использования ННТ и EMD при анализе сейсмической активности и землетрясений [8], анализе структуры сейсмических сигналов [9, 10], изучении сейсмоакустиче-ской эмиссии горных пород [11], подавлении помех в каротажных данных [12], обработке материалов аэрогравиметрических исследований [13].
В методе EMD предполагается, что анализируемый сигнал состоит из серии составляющих с различными частотами (intrinsic mode functions -IMF), обладающих следующими свойствами: 1) число максимумов и минимумов функции, а также и
количество пересечений нуля отличаются не более, чем на единицу; 2) среднее значение огибающих, построенных по локальным максимумам и локальным минимумам, близко к нулю. В отличие от гармоник, получаемых при представлении дискретного сигнала рядом Фурье, каждая IMF может иметь переменную амплитуду и частоту в разные моменты времени t (или в разных точках пространства x).
Для геофизических полей характерна нестационарность, т. е. естественное изменение их статистических характеристик в пространстве [14]. Природа этих изменений различна, и, в частности, может быть связана с фрактальными особенностями полей, обладающих самоподобной иерархически упорядоченной структурой [15]. EMD может широко применяться для анализа нестационарных данных, т. е. является адекватным геофизической практике, в отличие от большинства известных методов. Несмотря на еще недостаточно разработанное теоретическое обоснование, можно полностью согласиться с тем, что «метод эмпирических мод выступает в роли нового альтернативного (и также обладающего широкими возможностями) инструмента исследования структуры сложных сигналов» [16].
В данной статье предлагается модификация метода EMD, ориентированная на обработку геофизических данных. Отличительной особенностью ее является использование гармонических (потенциальных) функций при построении огибающих сигнала вместо традиционно применяющихся для этой цели кубических сплайнов.
Алгоритм EMD и его программная реализация
Возьмем значения /(x) геофизического поля (сигнал), полученные в k точках наблюдений на профиле с постоянным шагом Ax. Выделим все максимумы и минимумы сигнала и проведем через эти точки огибающие p (x) и q (x), соответственно. Теперь для всех k точек можно определить функцию средних значений <(x)=[p(x)+q(x)]/2. Разность V(x) между/(x) и <(x) будет являться первым приближением IMF 1 (рис. 1) - первой компонентой отсеивания (sifting).
/V
тически выполняемое число итераций составляет n=n(S). В результате отсеивания будет получена IMF 1. Для построения IMF 2 нужно будет организовать новый цикл (1), используя теперь в качестве исходных данных разность сигнала и уже имеющейся модовой функции: V2,i=/(x)-Vin Дальнейшее вычисление всех IMF (разложение: внешний итерационный цикл по i) проводится на основе последовательного исключения из сигнала всех предыдущих модовых функций:
V,+и = ~Wt, > 2 ^ * < m (3)
Для решения многих геофизических задач характерно наличие фоновой компоненты в анализируемом сигнале. Поэтому для останова внешнего цикла (3) используется критерий, предполагающий высокую гладкость остатка r(x)=ym+u(x): вычисления прекращаются, если функция r(x) имеет меньше двух максимумов и меньше двух минимумов. Таким образом, после проведения всех циклов просеивания (1) и разложения (3) исходный сигнал удается представить в виде суммы всех IMF и остаточной составляющей:
f ( x ) = ,п ( x) +r( x).
(4)
Рис. 1. Построение первого приближения IMF Vß 1 ~ сигнал fix); огибающие: 2 ~ р(х), 3 ~ q(x); функции: 4 ~ p(x), 5 - y(x)
Fig. 1. Construction of the first approximation to the IMF 1:1 is the f(x) signal; envelopes: 2 - p(x), 3 - q(x); functions: p(x), 5 - y(x)
Для дальнейшей записи введем нижние индексы i=1, 2,..., m, отвечающие вычисляемой IMF (а также функции (px)) и номеру j=1, 2,..., n ee приближения. Теперь охарактеризованное выше первое приближение к IMF 1 - это yu(x)=f(x)-pu(x). Процесс уточнения каждой IFM (отсеивание: внутренний итерационный цикл по j) выглядит следующим образом:
Wij+i(x) j (x) -p j (x), 1 < j < n. (1)
Критерием его завершения является выполнение условия:
X M j( x ) j-i( x )]2 < e
s >.,-i( x ) ]■ (2)
Сигнал ¡(х), согласно (4), оказывается разложенным по отвечающему исходным данным конечному адаптивному базису, не имеющему аналитического описания. Этот базис является полным, ортогональным и, по мнению Н. Хуанга, единственным [17].
Важнейшим элементом процесса эмпирической модовой декомпозиции является построение огибающих сигнала ^(х) и д(х). Для этой цели широко используются кубические сплайны, также апробировано применение В-сплайнов [18]. Однако методы сплайн-интерполяции нередко вносят заметные искажения в краевые части сигнала, для уменьшения которых в данной статье предлагается использование истокообразной аппроксимации при вычислении огибающих ^(х) и д(х). Алгоритм состоит в следующем: первоначально определяются координаты х точек локальных экстремумов анализируемой функцииДх) (или у(х)) и значения этих экстремумов. Предположим, что имеется 1Х значений амплитуды Нх максимумов и Ц значений амплитуды Н2 минимумов.
Можно рассчитать наибольшее расстояние хюах между соседними парами максимумов и аппроксимировать огибающую ^(х) линейной комбинацией гармонических функций:
p(x) = £b,z/[(x, - x)2 + z2],
(5)
где 6 - заданная погрешность или осуществление требуемого числа итераций п. Таким образом, фак-
где Ь коэффициенты, определяемые в процессе аппроксимации; г - константа; - абсциссы локальных максимумов. Правая часть выражения (5) представляет собой гравитационный эффект совокупности бесконечно длинных горизонтальных стержней, размещенных на глубине г под точками
i=i
t=i
xt, i=l,2,^,i1, обладающих линейными массами, пропорциональными bt [19]. Значения коэффициентов bt определяются путем решения системы линейных алгебраических уравнений (СЛАУ) методом Зейделя:
Ab=c, (6)
где A - матрица коэффициентов вида z/[(xt-x)2+z2] размером t1 строк, t1 столбцов; b={bt} - вектор неизвестных параметров; c={hj - вектор максимальных значений функции f(x) или y/(x). Выбор константы 0,5xmax<z<2xmax обеспечивает устойчивость вычислений и высокую точность аппроксимации [20, 21]. Сходимость метода обусловлена тем, что матрица A является симметричной, положительно определенной и наделена свойством диагонального преобладания. Качество аппроксимации контролируется в метрике Чебышева:
max \hl( xt) - p( xt) | <s0, (7)
где s0 - достаточно малая положительная величина.
Построение огибающей по минимумам q(x) осуществляется аналогично, в соответствии с формулами (5)-(7), при замене в них p (x) на q (x), t1 на t2, h1 на h2. Полученные аналитические выражения (5) для двух огибающих используются при вычислении функции средних значений ^(x) во всех к точках оси x. Высокую точность построения огибающих p(x) и q(x), а также незначительные искажения в краевых частях профиля иллюстрирует рис. 1.
Вышеописанный алгоритм эмпирической модо-вой декомпозиции реализован авторами в программе RIMF, созданной с использованием системы визуального объектно-ориентированного программирования DELPHI 8.0. В программе используются следующие параметры: максимальное количество эмпирических мод т<15; показатель качества аппроксимации s0<104L, где L - размах аппроксимируемой функции; пользователем задаются максимальное число циклов отсеивания n (по умолчанию n<7 и погрешность (3) приближенного построения IMF (по умолчанию 5<0,01. Как правило, для достижения условия (7) требуется не более 20-30 итераций при решении СЛАУ (6). Затраты времени на декомпозицию данных аэромагнитной съемки, заданных в к=601 точках, при i=7 и суммарном количестве циклов просеивания Zn=43, на компьютере с процессором Intel (R) Core (TM) i7-4770K с тактовой частотой 3,5 ГГц составляют ~10 с.
Примеры применения метода
Рассмотрим фильтрационные возможности алгоритма EMD на примере данных электропрофилирования методом срединного градиента, с питающей линией AB=400 м и приемной линией MN=10 м, проведенного с целью решения инженерно-геологических задач в пределах Верхнекамского месторождения калийно-магниевых солей (Пермский край). Профиль наблюдений включает 138 точек, шаг между точками 10 м. Полученные в результате полевых измерений значения кажущегося электрического сопротивления рк, лежа-
щие в диапазоне от 0,7 до 129 Омм (рис. 2, а), в результате декомпозиции были разложены на шесть IMF и остаток r(x). С целью подавления помех использовалось усреднение в скользящем окне размером 7 точек (рис. 2, б) и синтез полезного сигнала (рис. 2, в), отождествляющегося в данном случае с суммой IMF с номерами 3, 4, 5, 6 и остаточной компоненты r(x). Как очевидно, в последнем случае результативный график является более гладким, при этом не происходит потери информации на участках длиной 30 м в начале и конце профиля.
О 0.4 0.8 1.2 X. км
Рис. 2. Графики кажущегося электрического сопротивления: а) результаты измерений; сглаженные значения, полученные с помощью: б) осреднения в скользящем окне; в) суммирования нескольких IMF и остатка r(x). Примечание: черный пунктир - остаток r(x)
Fig. 2. Graphs of apparent electrical resistance: а) the results of the measurements; smoothed values are obtained by: б) averaging in a sliding window; в) using summation of multiple IMFs and the residue r(x). Note: the remainder r(x) is marked with the black dashed line
Метод эмпирической модовой декомпозиции можно применять для приближенной оценки уровня аппаратурных и геологических помех в результатах геофизических съемок. Приведем пример подобной оценки по результатам высокоточной гравиметрической съемки масштаба 1:25000, выполненной в Западном Саяне гравиметрами Au-tograv CG-5 (Scint^, Канада) с привязкой гравиметрических пунктов при помощи спутниковой системы GPS Trimble-4700 и электронных тахео-
Рис. 3. Аномалии силы тяжести и их региональная компонента (а); выделенные методом EMD составляющие, отвечающие: погрешности гравиметрической съемки (б), помехам геологической природы (в), изучаемым геоплотностным неодно-родностям (г)
Fig. 3. Gravity anomalies and their regional component (а); components allocated by the EMD method meet: the error of the gravity survey (б), the nature of geological objects (в), the geo-density inhomogeneities (г)
метров с целью поисков медно-молибденового ору-денения. Интерполяционный гравиметрический профиль, пересекающий участок работ в широтном направлении, включает в себя 200 точек измерений. Расстояние между точками 100 м. Перепад значений аномального гравитационного поля Ag составляет более 7 мГал, его значения осложнены региональным фоном, для исключения которого использовалась аппроксимация полиномом 2-й степени: Ag[1(i=27,535+0,069x+0,014x2 (рис. 3, а). Локальная составляющая поля содержит 5 эмпирических мод. IMF 1 - это наиболее высокочастотная составляющая локальной компоненты поля Ag^Ag-Agp,,,,, обусловленная погрешностью съемки (рис. 3, б). Ее статистическое распределение близко к нормальному (рис. 4), математическое ожидание M=0,002, среднеквадратическое отклонение (СКО) составляет ±0,067 мГал. Эта величина близка к среднеквадратической погрешности определения аномалии Буге, которая составляет ±0,092 мГал (с учетом погрешностей наблюдений, планово-высотной привязки и вычисления поправок за влияние рельефа местности). Сумма двух следующих IMF: предположительно отра-
жает влияние помех геологического происхождения (рис. 3, в). СКО этой компоненты ±0,091 мГал. Полезный сигнал, представленный на рис. 3, г, синтезирован путем суммирования остальных компонент разложения локальной составляющей поля щ+щ+r.
Рис. 4. Статистическое распределение значений IMF 1 и его аппроксимация законом Гаусса (черный пунктир)
Fig. 4. Statistical distribution of IMF 1 values and its approximation by the Gaussian (black dotted line)
При моделировании геологического строения по ряду профилей с использованием монтажного метода решения обратной задачи гравиразведки [21] среднеквадратическое расхождение наблюденного и модельного полей составляло 0,2-0,25 мГал, что хорошо согласуется с представленными выше оценками помех.
Скрытые особенности взаимосвязи геофизических полей и аномалиеобразующих геологических объектов можно пытаться выявить на основе применения метода EMD. В определенных диапазонах пространственных частот эти особенности могут проявляться достаточно отчетливо, поэтому целесообразно проводить сопоставление IMF, характеризующих аномалии и их источники. Рассмотрим один пример такого сопоставления, выполненного с целью приближенного учета влияния рельефа земной поверхности на результаты крупномасштабной аэромагнитной съемки, выполненной в Норильском районе НФ ВСЕГЕИ в 2012-2013 гг.
Широко развитая на территории съемки туфо-лавовая толща, максимальная мощность которой превышает 3000 м, является главной помехой при поисках рудоносных интрузий геофизическими методами [23]. По направлению (знаку) естественной остаточной намагниченности In туфолавовая толща разделяется на три неравномерные части. Нижняя отвечает только породам ивакинской свиты. Верхняя включает два верхних потока базальтов самоедской свиты, намагниченных отрицательно. Остальная часть разреза намагничена положительно. Суммарная намагниченность 1эф ту-фолавовой толщи Норильского района равна 3з3х10-2А/м при следующих значениях средневзвешенной суммарной намагниченности отдельных свит, А/м: ивакинская - 40х10-2, сывермин-ская - 43х10-2, гудчихинская вместе с хаканча-нской - 107х10-2, туклонская - 78х10-2, надеждин-ская - 192х10-2, моронговская - 412х10-2, мокула-евская - 207х10-2, хараелахская - 616х10-2, кум-гинская - 629х10-2, самоедская - 527х10-2. Влияние резко расчлененного рельефа местности, сложенного интенсивно намагниченными горными породами, выражается в пространственной корреляции повышений амплитуды аномального магнитного поля с увеличением высотных отметок. Однако, в силу сложной физико-геологической обстановки, зависимость магнитного поля от высоты носит нелинейный характер и для его редуцирования используются достаточно сложные технологии, аналогичные вычислению поправок за влияние рельефа в гравиразведке. В разные годы Г.Г. Ремпелем (СНИИГГиМС), В.Н. Юровских (СЕГФЭ), П.В. Кирплюком (НФ ВСЕГЕИ), А.С. Долгалем (Таймыргеолком) были разработаны оригинальные методы определения магнитных аномалий, обусловленных рельефом дневной поверхности, приближенно учитывающие латеральную изменчивость намагниченности трапповых образований.
В приведенном примере значения аномального геомагнитного поля AT заданы по широтному профилю длиной 60 км, пересекающему Талнахское месторождение медно-никелевых руд. Амплитуда поля меняется от -627 до 512 нТл, шаг между точками 100 м (рис. 5).
В результате модовой декомпозиции для магнитного поля над туфолавовой толщей было полу-
чено 7 IMF, для высотных отметок рельефа -4 IMF. Характеристика этих функций представлена в виде диаграмм Тьюки («ящик с усами»), приведенных на рис. 6. Взаимосвязь между модальными компонентами поля и рельефа отражает таблица парных коэффициентов линейной корреляции (таблица).
Рис. 5. Аномальное магнитное поле (а) и рельеф земной поверхности (б): 1 - осадочные породы; 2 - туфолавовая толща; 3 - Талнахское месторождение
Fig. 5. Anomalous magnetic field (a) and topography (b): 1 -are the sediments; 2 - strata of tuffs and lavas; 3 is the Talnakh Deposit
Более 50 % коэффициентов корреляции K не являются значимыми. Значимость оценивалась с помощью í-критерия Стьюдента для доверительной вероятности 99,9 % (v=339). В целом между высотами H и магнитным полем AT отмечается слабая корреляционная связь (K=0,295), однако взаимосвязь между IMF этих параметров является более тесной (значения K достигают 0,976). Это позволяет построить регрессионные уравнения вида \\(AT)=A+B\(H), 1<i<7, 1<j<4 и рассчитать с их помощью 7 разночастотных составляющих магнитного поля, связанных с изменением высот земной поверхности. Сумма этих составляющих по своему физическому смыслу будет отвечать влиянию неоднородно намагниченного рельефа туфолавовой толщи, т. е. являться поправкой за влияние рельефа STp. Полученные значения STp достаточно хорошо согласуются с независимыми результатами вычисления поправок за влияние рельефа местности, полученными путем решения прямой задачи магниторазведки (рис. 7). При этом использовалась детальная цифровая модель высот рельефа с переменной по латерали намагниченностью J=J(x,y), радиус учета влияния рельефа составлял 20 км [24]. Следует отметить, что оба способа расчета поправок STp являются приближенными, но предложенный авторами способ, базирующийся на EMD, требует значительно меньших трудозатрат.
IMF 1 IMF 2 IMF 3 IMF 4 IMF 5 IMF 6 IMF7
Рис. 6. Диаграммы Тьюки для магнитного поля (а) и высот рельефа местности (б) Fig. 6. Tukey diagrams for the magnetic field (а) and heights of the terrain (б)
Таблица. Парные коэффициенты линейной корреляции K между магнитным полем AT, высотами рельефа H и компонентами их разложения ц, r
Table. Paired coefficients of linear correlation K between the magnetic field AT, the height H and the components of decomposition
Параметры/Parameters AT Щ Щ Щз Щ4 Щ5 Щ6 Щ r
H 0,295 0,078 0,092 0,171 0,314 0,433 0,450 0,209 -0,076
Щ 0,224 0,644 0,141 0,0193 -0,030 -0,095 -0,104 -0,129 -0,043
Щ 0,338 0,082 0,614 0,235 -0,117 0,011 0,028 0,012 0,075
Щ 0,317 -0,029 0,110 0,303 0,083 0,354 0,304 -0,014 0,184
Щ4 0,546 -0,105 -0,076 0,168 0,727 0,976 0,947 0,674 0,167
r -0,092 0,052 -0,053 -0,026 0,093 0,009 0,063 -0,003 -0,268
Примечание: жирным шрифтом выделены компоненты магнитного поля AT, для которых составлены уравнения линейной регрессии.
Note: the components of the magnetic field AT with the regression equations are in bold.
Рис. 7. Магнитное поле STp над туфолавовой толщей, обусловленное влиянием рельефа местности, определенное с использованием метода EMD (красная линия) и решения прямой задачи магниторазведки (синяя пунктирная линия)
Fig. 7. Magnetic field 8Tf above the strata of tuffs and of lavas, due to the influence of the terrain, determined using the EMD method (red line) and solution of the direct problem of magnetometry (blue dotted line)
Заключение
Представленные в статье теоретические и практические материалы позволяют сделать вывод о целесообразности применения метода эмпи-
рическом модовой декомпозиции при анализе геофизических данных. Использование адаптивного базиса разложения сигнала в ряде случаев позволяет выделить его физически значимые
компоненты, отвечающие влиянию отдельных факторов. В частности - наиболее высокочастотную компоненту гравитационного поля (IMF 1) можно сопоставить с погрешностью определения аномалий Буге. Выделение информативных составляющих осложненных помехой геофизических полей осуществляется путем суммирования нескольких модовых функций щ с номерами 2<i<m, где m - номер последней IMF, полученной в процессе разложения сигнала. В данной реализации алгоритма EMD остаток r отвечает региональному фону (например, влиянию глубинных источников поля).
СПИСОК ЛИТЕРАТУРЫ
1. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis / N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, H.H. Shih, Q. Zheng, N.-C. Yen, C.C. Tung, H.H. Liu // Proc. R. Soc. Lond. A. - 1998. - Т. 454. -P. 903-995.
2. Huang N.E., Samuel S.S.P. Hilbert-Huang Transform and its applications. - Singapore: World Scientific Publishing Co, 2005. - 323 p.
3. Ястребов И.П. О свойствах и применениях преобразования Гильберта-Хуанга // Проектирование и технология электронных средств. - 2016. - № 1. - С. 26-33.
4. Нежевенко Е.С., Феоктистов А.С. Преобразование Гильберта-Хуанга двумерных изображений и использование его для выделения текстурных признаков // Региональные проблемы дистанционного зондирования Земли: Материалы международной научной конференции. - Красноярск, Россия, 23-26 сентября 2014. - Красноярск: Сиб. федер. ун-т, 2014. - С. 200-203.
5. Сафиуллин Н.Т. Разработка методики анализа временных рядов с помощью преобразования Хуанга-Гильберта: дис. ... канд. техн. наук. - Новосибирск, 2015. - 193 с.
6. Kahaner D., Moler C., Nash S. Numerical Methods and Software Cleve Moler and Stephen Nash. - Englewood Clis: Prentice-Hall, 1989. - 495 p.
7. Dremin I.M., Ivanov O.V., Nechitailo V.A. Wavelets and their uses // Physics-Uspekhi. - 2001. - V. 44. - № 5. - P. 447-478.
8. Huang N.E., Wu Z. A review on Hilbert-Huang transform: Method and its applications to geophysical studies // Rev. Geo-phys. - 2008. - 46. - RG2006. DOI: 10.1029/2007RG000228.
9. Павлов А.Н., Филатова А.Е. Метод эмпирических мод и вей-влет фильтрация: применение в задачах геофизики // Известия вузов. Прикладная нелинейная динамика. - 2011. -Т. 19. - № 1. - С. 3-13.
10. Долгих Л.А. Выделение информативной составляющей сейсмо-акустического сигнала: возможные подходы // Проблемы автоматизации и управления в технических системах: Тр. между-нар. науч. техн. конф. - Пенза: ИИЦ ПГУ, 2008. - С. 322-323.
11. Иванов Д.Б., Бабуркина М.А. Анализ изменения частотного состава сигналов естественной и вызванной сейсмоакустиче-ской эмиссии // Двенадцатая уральская молодежная научная школа по геофизике: Сборник науч. материалов. - Пермь: ГИ УрО РАН, 2011. - С. 86-90.
12. Давыдов В.А., Давыдов А.В. Управление эмпирической модо-вой декомпозицией сигналов при анализе и обработке геофизических данных // Каротажник. - 2010. - № 5. - С. 98-114.
Представляет интерес применение предварительного выделения модовых функций в корреляционных методах интерпретации геофизических данных. Это подтверждают удовлетворительные результаты определения поправок за влияние рельефа на магнитное поле STp, полученные в сложных физико-геологических условиях Норильского района. Значительными перспективами обладают методы EMD и HHT при анализе материалов геофизического мониторинга месторождений полезных ископаемых, находящихся в стадии разработки.
Работа выполнена при поддержке грантов РФФИ № 16-45-590046, № 15-05-01823.
13. Hassan H.H., Peirce J.W. Empirical Mode Decomposition (EMD) of potential field data: airborne gravity data as an example // Recorder. - January 2008. - V. 33. - № 01. - P. 25-30. URL: http://csegrecorder.com/articles/view/empirical-mode-decom-position-emd-of-potential-field-data/ (дата обращения 27.10.2016).
14. Никитин А.А, Петров А.В. Теоретические основы обработки геофизической информации. - М.: РГГУ, 2008. - 112 с.
15. Блох Ю.И. Проблема адекватности интерпретационных моделей в гравиразведке и магниторазведке // Геофизический вестник. - 2004. - №6. - С. 10-15.
16. Ульянова Ю.Е., Бабенко Р.Г., Чернов А.В. Частотно-временные преобразования, используемые в цифровой обработке сигналов // Глобальная ядерная безопасность. - 2015. - № 3 (16). - С. 36-42.
17. Давыдов В.А., Давыдов А.В. Очистка геофизических данных от шумов с использованием преобразования Гильберта-Хуанга // Электронное научное издание «Актуальные инновационные исследования: наука и практика». - 2010. - № 1. URL: http://www.actualresearch.ru/nn/2010_1/Article/geo/davy-dov.htm (дата обращения 17.10.2016).
18. Chen Q., Huang N., Riemenschneider S. A B-spline approach for empirical mode decomposition // Adv. Comp. Math. - 2004. -№24. - P. 171-195.
19. Гравиразведка: справочник геофизика / под ред. Е.А. Мудре-цовой, К.Е. Веселова. - М.: Недра, 1990. - 607 с.
20. Strakhov V.N., Stepanova I.E. The S-Approximation Method and its Application to Gravity Problems // Izvestiya, Physics of the Solid Earth. - 2002. - V. 38. - № 2. - P. 91-107.
21. Effective Algorithms for Sourcewise Approximation of Geopo-tential Fields / P.I. Balk, A.S. Dolgal, A.V. Pugin, A.V. Michu-rin, A.A. Simanov, A.F. Sharkhimullin // Izvestiya, Physics of the Solid Earth. - 2016. - V. 52. - № 6. - P. 896-911.
22. Balk P.I., Dolgal A.S. Three-dimensional assembly technologies for the interpretation of gravimetric data // Doklady Earth Sciences. - 2009. - V. 427. - № 2. - P. 971-974.
23. Dolgal A.S., Chekhovich K.M. Complex interpretation of geopoten-tial fields on searches or copper-nickel-platinum metallization (Norilsk region) // Geologiya i geofizika. - 1998. - V. 39. - № 11. -P. 1615-1625.
24. Долгаль А.С. Магниторазведка: компьютерные технологии учета влияния рельефа местности. - Пермь: Перм. гос. нац. ис-след. ун-т, 2014. - 92 с.
Поступила 10.01.2017 г.
Информация об авторах
Долгаль А.С., доктор физико-математических наук, главный научный сотрудник лаборатории геопотенциальных полей Учреждения Российской академии наук «Горный институт Уральского отделения РАН».
Христенко Л.А., кандидат геолого-минералогических наук, научный сотрудник лаборатории наземной и подземной электрометрии Учреждения Российской академии наук «Горный институт Уральского отделения РАН».
UDC 550.83.016
APPLICATION OF EMPIRICAL MODE DECOMPOSITION METHOD IN PROCESSING GEOPHYSICAL DATA
Alexander S. Dolgal1,
Liudmila A. Khristenko1,
1 Mining Institute of the Ural Branch Russian Academy of Sciences, 78-a, Sibirskaya street, Perm, 614007, Russia.
The relevance of the research is caused by the feasibility of using advanced mathematical methods in processing results of the geophysical surveys.
The aim of the research is to improve the effectiveness of the methods of applied Geophysics through adaptive extraction of informative components of the physical fields using the Empirical Mode Decomposition method (EMD).
The method of research: decomposition of the relevant geophysical data. The iterative computational process allows allocating different frequency orthogonal signal components, which are called empirical mode functions (iMFs). The method is designed to represent non-stationary signals in the form of a series of signals with different frequency The authors propose the original algorithm in which the sourcewise approximation is used for constructing functions, enveloping the extrema of the analyzed signal, and the residual component of decomposition can identify the background component of the field.
The results. The EMD algorithm was implemented and tested on practical materials for processing the data of gravity, magnetic and electrical prospecting. The authors defined the possibility of applying the method to separate the informative component of the results of electric methods of horizontal profiling (Perm Krai); to evaluate the quality of high-precision gravity survey (Western Sayan); to determine the component of the magnetic field caused by the impact of sharply dissected terrain of the earth surface in the development of the effusive traps (Norilsk region). The last case uses a series of regression dependencies between the various frequency components of the magnetic field and heights. The authors made the conclusion on appropriateness of using the EMD technique in applied Geophysics. The prospects of its application the authors see in analysis of geophysical monitoring material of mineral deposits development.
Key words:
Signal, Empirical Mode Decomposition, algorithm, approximation, field component, gravity survey, magnetic survey, electrical prospecting.
The research was supported by the grants of RFBR no. 16-45-590046, no. 15-05-01823.
REFERENCES
1. Huang N.E., Shen Z., Long S.R., Wu M.C., Shih H.H., Zheng Q., Yen N.-C., Tung C.C., Liu H.H. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. Lond. A., 1998, vol. 454, pp. 903-995.
2. Huang N.E., Samuel S.S.P. Hilbert-Huang Transform and its applications. Singapore, World Scientific Publishing Co, 2005. 323 p.
3. Yastrebov I.P. On properties and applications of Hilber-Huang transform. Proektirovanie i tekhnologiya elektronnykh sredstv, 2016, no. 1, pp. 26-33. In Rus.
4. Nezhevenko E.S., Feoktistov A.S. Preobrazovanie Gilberta-Huan-ga dvumernykh izobrazheniy i ispolzovanie ego dlya vydeleniya teksturnykh priznakov [The Hilbert-Huang transform of two-dimensional images and its use to highlight texture features]. Regio-nalnye problemy distantsionnogo zondirovaniya Zemli. Materialy mezhdunarodnoy nauchnoy konferentsii [Materials of the international scientific conference. Regional problems of Earth remote sensing]. Krasnoyarsk, Russia, 23-26 September 2014. Krasnoyarsk, Siberian Federal University Press, 2014. pp. 200-203.
5. Safiullin N.T. Razrabotka metodiki analiza vremennykh ryadov s pomoshchyu preobrazovaniya Huanga-Gilberta. Dis. kand. nauk [Development of methods of time series analysis by the Huang-Hil-bert converting. Cand. Diss.]. Novosibirsk, 2015. 193 p.
6. Kahaner D., Moler C., Nash S. Numerical Methods and Software Cleve Moler and Stephen Nash. Englewood Clis, Prentice-Hall, 1989. 495 p.
7. Dremin I.M., Ivanov O.V., Nechitailo V.A. Wavelets and their uses. Physics, Uspekhi, 2001, vol. 44, no 5, pp. 447-478.
8. Huang N.E., Wu Z. A review on Hilbert-Huang transform: Method and its applications to geophysical studies. Rev. Geophys., 2008, 46, RG2006. DOI: 10.1029/2007RG000228.
9. Pavlov A.N., Filatova A.E. Empirical mod and wavelet filtering: application to Geophysics problems. Izvestiya vuzov. Prikladnaya nelineynaya dinamika, 2011, vol. 19, no. 1, pp. 3-13. In Rus.
10. Dolgikh L.A. Vydelenie informativnoy sostavlyayushchey seys-moakusticheskogo signala: vozmozhnye podkhody [Allocation of the informative component of the acoustic signal: possible approaches]. Problemy avtomatizatsii i upravleniya v tekhnicheskikh si-stemakh: Trudy mezhdunarodnoy nauchnoy tekhnicheskoy konfe-rentsii [Proceedings of the international scientific technical conference. Problems of automation and control in technical systems]. Penza, IIC PGU, 2008, pp. 322-323.
11. Ivanov D.B., Baburkina M.A. Analiz izmeneniya chastotnogo sostava signalov estestvennoy i vyzvannoy seysmoakustiches-koy emissii [Analysis of changes of frequency composition of signals of natural and induced seismoacoustic emissions]. Dven-adtsataya uralskaya molodezhnaya nauchnaya shkola po geofi-zike: Sbornik nauchnykh Materialov [Twelfth Ural youth scientific school on Geophysics: Proc.]. Perm, Mining Institute of the Ural Branch of the Russian Academy of Sciences, 2011. pp. 86-90.
12. Davydov V.A., Davydov A.V. Management of Empirical Mode Decomposition of signals in the analysis and processing of geophysical data. Karotazhnik, 2010, no. 5, pp. 98-114. In Rus.
13. Hassan H.H., Peirce J.W. Empirical Mode Decomposition (EMD) of potential field data: airborne gravity data as an example. Recorder, January 2008, vol. 33, no. 01, pp. 25-30. Available at: http://csegrecorder.com/articles/view/empirical-mode-decom-position-emd-of-potential-field-data/ (accessed 27 October 2016).
14. Nikitin A.A, Petrov A.V. Teoreticheskie osnovy obrabotki geofi-zicheskoy informatsii [Theoretical bases of processing geophysical data]. Moscow, RGGU Press, 2008. 112 p.
15. Blokh Yu.I. Problema adekvatnosti interpretatsionnykh modeley v gravirazvedke i magnitorazvedke [The problem of adequacy of the interpretive models in gravity and magnetic]. Geofizicheskiy vestnik, 2004, no. 6, pp. 10-15.
16. Ulyanova Yu.E., Babenko R.G., Chernov A.V. The time frequency conversion used in digital signal processing. Globalnaya yader-naya bezopasnost, 2015, no. 3 (16), pp. 36-42. In Rus.
17. Davydov V.A., Davydov A.V. Cleaning geophysical data from noise using the Hilbert-Huang transform. Elektronnoe nauchnoe iz-danie «Aktualnye innovatsionnye issledovaniya: nauka i praktika», 2010, no. 1. Available at: http://www.actualresearch.ru (accessed 17.10.2016). In Rus.
18. Chen Q., Huang N., Riemenschneider S. A B-spline approach for empirical mode decomposition. Adv. Comp. Math., 2004, no. 24, pp.171-195.
19. Gravirazvedka: spravochnik geofizika [Gravity: geophysics guide]. Eds. E.A. Mudretsova, K.E. Veselov. Moscow, Nedra Publ., 1990. 607 p.
20. Strakhov V.N., Stepanova I.E. The S-Approximation Method and its Application to Gravity Problems. Izvestiya, Physics of the Solid Earth, 2002, vol. 38, no. 2, pp. 91-107.
21. Balk P.I., Dolgal A.S., Pugin A.V., Michurin A.V., Sima-nov A.A., Sharkhimullin A.F. Effective Algorithms for Sour-cewise Approximation of Geopotential Fields. Izvestiya, Physics of the Solid Earth, 2016, vol. 52, no. 6, pp. 896-911.
22. Balk P.I., Dolgal A.S. Three-dimensional assembly technologies for the interpretation of gravimetric data. Doklady Earth Sciences, 2009, vol. 427, no. 2, pp. 971-974.
23. Dolgal A.S., Chekhovich K.M. Complex interpretation of geopotential fields on searches or copper-nickel-platinum metallization (Norilsk region). Geologiya i geofizika, 1998, vol. 39, no. 11, pp. 1615-1625.
24. Dolgal A.S. Magnitorazvedka: kompyuternye tekhnologii ucheta vliyaniya relefa mestnosti [Magnetic: computer technology of considering the influence of the terrain]. Perm, Perm State National Research University Press, 2014. 92 p.
Received: 10 January 2017.
Information about the authors
Alexander S. Dolgal, Dr. Sc., chief research worker, Mining Institute of the Ural Branch Russian Academy of Sciences.
Liudmila A. Khristenko, Cand. Sc., researcher, Mining Institute of the Ural Branch Russian Academy of Sciences.
