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

УДК 550.37 Д.В. Епишкин1

РАЗВИТИЕ МЕТОДОВ ОБРАБОТКИ ДАННЫХ МАГНИТОТЕЛЛУРИЧЕСКОГО ЗОНДИРОВАНИЯ

Разработан алгоритм (система) обработки магнитотеллурических (МТ) данных, обладающий высокой устойчивостью по отношению к интенсивным электромагнитным шумам, присутствующим в записях магнитотеллурического поля. Особенность алгоритма — применение специальных схем мультипликативного взвешивания данных при оценивании разных передаточных функций. В алгоритме использованы разные методы борьбы со смещением спектральных оценок, включая робастную М-оценку, метод «складного ножа», использование удаленной базы, ввод поправок за смещение спектральных оценок. Предложенный алгоритм демонстрирует высокую эффективность при обработке МТ-данных с низким соотношением сигнал/шум.

Ключевые слова: обработка магнитотеллурических данных, оценка передаточных функций, робастное оценивание.

A magnetotelluric data processing code has been developed, which demonstrates high robustness to intense electromagnetic noise occurring in measured MT data. Key features of the code are specific approach for estimating different transfer functions and capability to utilize all four channels acquired at remote reference station. The code utilizes various techniques to reduce estimate errors, including robust Huber estimator, jackknife approach, improved remote reference technique and compensating for overestimation of power spectra. The proposed code has shown high efficiency in processing of low signal-to-noise data.

Key words: magnetotelluric data processing, response function estimation, robust estimator.

Введение. При проведении магнитотеллурических (МТ) исследований в производственных масштабах в большинстве российских организаций используются аппаратура и программное обеспечение зарубежных производителей. Если ориентироваться на развитие технологической базы, в том числе на разработку новых образцов отечественной аппаратуры для магнитотеллурического зондирования (МТЗ), то необходимо создать и соответствующее программное обеспечение. При этом важная задача — повышение эффективности процедур обработки по сравнению с широко применяемыми на практике в настоящее время, особенно в условиях низких значений отношения сигнал/шум.

Указанные предпосылки определяют актуальность создания новых вычислительных алгоритмов и программных средств для обработки МТ-данных, т.е. для оценки магнитотеллурических передаточных операторов. Предлагается алгоритм обработки МТ-данных, учитывающий особенности различных передаточных операторов; опробованы подходы к уменьшению смещения оценок.

Задачи, принципы и проблемы обработки МТ-данных. Под обработкой данных МТЗ понимают нахождение передаточных функций (тензор импеданса, матрица Визе—Паркинсона, теллурический и горизонтальный магнитный тензор) по записям

компонент магнитотеллурического поля. При этом основной этап обработки — разделение записей МТ-поля на спектральные составляющие, по которым затем находятся компоненты искомых передаточных функций для заданного набора частот [Семенов, 1985; Larsen, 1989].

Базовый подход, ограничивающийся анализом записей в точке зондирования, подробно изложен, например, в [Семенов, 1985; Sims, 1971]. В его основе лежит поиск такой передаточной функции, например, тензора импеданса

Z —

7

v

ху

"УУ J

Ex Zxx ' Hx + ZxyHy, Ey = Zyx ' Hx + ZyyHy,

которая удовлетворяла бы (в среднеквадратичном или каком-либо ином смысле) всему множеству реализаций наблюденного МТ-поля:

(1а) (1б)

где Ех, Еу — спектры горизонтальных компонент электрического поля; Нх, Ну — спектры горизонтальных компонент магнитного поля соответственно. Все величины в уравнениях (1а, 1б) — функции частоты. Уравнения (1) принято решать, разбивая запись на сегменты (окна), выполняя для них преобразование Фурье и затем минимизируя

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

Рис. 1. Удаление «пиков» в записях электрического поля. Черное — пики, серое — исправленный временной ряд

невязку между левой и правой частями системы (1), например, с помощью метода наименьших квадратов (МНК). Недостаток стандартного МНК заключается в существенной зависимости результата от наличия и характеристик шума, всегда присутствующего в измеряемых компонентах поля. Предложен ряд усовершенствований этого подхода, включающих различные статистические робастные процедуры [Chave, 1987; Egbert, 1986], а также метод удаленной базовой точки [Gamble, 1979]. Тем не менее в условиях низкого отношения сигнал/шум результаты их применения остаются подвержены искажениям (смещению оценок), часто очень значительным.

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

Предварительная обработка временных рядов. На первом этапе выполняется ряд преобразований исходных записей поля, направленных на удаление двух специфических видов шума, которые часто присутствуют в записях компонент электрического поля — интенсивных «пиков», значимо влияющих на спектр сигнала (рис. 1), а также «ступеней» (рис. 2).

Для поиска «пиков» выбирается окно (фрагмент временного ряда) некоторой ширины, и для него находятся медиана и робастная оценка дисперсии [Huber, 1981]. Если отсчет сигнала в центре окна отличается от медианы больше чем на заданное число дисперсий, то отсчет считается аномальным. Затем путем последовательной проверки соседних отсчетов определяется ширина аномальной зоны, попадающие в нее отсчеты сигнала заменяются интерполированными значениями.

Рис. 2. Коррекция «ступеней» в записях электрического поля. Серое — временной ряд до коррекции, черное — ряд после коррекции

Поиск «ступеней» происходит на основании двух критериев — большого (относительно других точек) значения модуля производной и значительного увеличения дисперсии сигнала в окрестности ступени. Алгоритм поиска аномалий в производной аналогичен алгоритму поиска пиков в сигнале. Начало ступени находится по максимуму производной (осредненной в окне из нескольких точек) и по максимуму дисперсии.

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

Расчет спектров и весовых функций. Расчет спектров выполняется в окне фиксированной ширины. Окно движется по временному ряду с заданным шагом, и для каждого его положения выполняется оконное преобразование Фурье. В результате применения такой схемы можно получить некоторое число оценок спектра каждой компоненты поля для набора частот. Затем определяется качество каждой спектральной оценки, а в соответствии с этим задается ее вес, т.е. число в диапазоне от 0 до 1.

Алгоритм использует две робастные процедуры для расчета весов:

— отбраковка по когерентности и невязке на основе метода «складного ножа» [Jones, 1989]. Выполняется поиск и удаление спектральных оценок, присутствие которых существенно ухудшает результат (уменьшает когерентность, увеличивает доверительные интервалы);

— робастная М-оценка (оценка Хьюбера [Huber, 1981]), которая в применении к обработке МТ-данных обсуждалась в работах [Chave, 1987; Egbert, 1986; Egbert, 1996]. Помимо оригиналь-

I.

ной весовой функции Хьюбера в описываемом алгоритме используется также функция Томпсона [Chave, 1987].

Отбраковка по критерию когерентности выполняется в первую очередь и является вспомогательной процедурой, способствующей более эффективному использованию М-оценки на последующем этапе. Отбраковка по когерентности — полезный инструмент даже в случае, когда соотношение сигнал/шум достаточно мало. В то же время робастная М-оценка оказывается эффективной в условиях слабого сигнала (так называемые мертвые диапазоны), а также в области низких частот, где выборки имеют небольшой размер. Отбраковка оценок с применением метода «складного ножа» дает лучшие результаты на высоких частотах, для которых размер выборок больше.

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

Например, применительно к тензору импеданса используются две оценки, отвечающие решениям методом взвешенных наименьших квадратов для двух формулировок задачи минимизации невязки уравнений (1) [Sims, 1971]:

-'ух

7

у ух

ху

"УУ J

-ху

-'УУ )

кЕуНх

ЕХЕХ

уЕуЕх

EyHyJ

ЕхЕу

н м

rV1

EyEyJ

уНуНх

НХЕХ

HyHyJ

уНуЕх

НхЕу

HyE*j

, (2a)

(2б)

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

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

м

(3)

1=1

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

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

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

Разработанный алгоритм можно применять к расчету всех основных передаточных функций: тензора импеданса, вектора Визе (типпера), теллурического и магнитного тензоров. В каждом случае весовые функции подбираются так, чтобы обеспечить максимальное подавление шума в соответствующих компонентах МТ-поля. Для каждой передаточной функции существуют особенности обработки.

Рассмотрим алгоритм оценивания теллурического тензора. В предлагаемом алгоритме задействованы временные ряды для 8 компонент МТ-поля. Помимо электрических компонент поля в базовой и рядовой точке используются магнитные компоненты в этих точках.

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

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

2) вычисление теллурического тензора для случая, когда рядовая точка полагается удаленной базой, а удаленная база — рядовой точкой. В этом случае действует обратное предположение — весь шум присутствует на удаленной базе;

3) вычисление тензора импеданса в рядовой точке. Веса, найденные в ходе этого вычисления, уменьшают влияние шума в электрических каналах рядовой точки;

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

При проведении опытных работ (Казахстан, 2014) установлено, что включение в обработку магнитных каналов существенно улучшает качество оценивания. На рис. 3 показан детерминант теллурического тензора как функция частоты. Приведено сравнение двух типов обработки — базового (при взвешивании используются процедуры, описанные в пунктах 1 и 2) и модифицированного (с добавлением пунктов 3 и 4). На каждой частоте получено 20 решений (весь временной ряд предварительно разбит на 20 сегментов, и для каждого была проведена независимая обработка в обоих вариантах). Очевидно, что дисперсия решений в случае задействования магнитных каналов меньше, особенно это заметно в интервале времени от 10 до 100 с.

Arg(rDET)

б

Рис. 3. Детерминант теллурического тензора как функция частоты, амплитуда и фаза. Сравнение двух типов обработки: а — стандартная обработка — серые точки; б — обработка с использованием магнитных каналов, черные точки

Описанный алгоритм оценки Т разработан для решения конкретной задачи — получения информации о поляризуемости пород [2огт, 2015]. Необходимо было получать фазу теллурического тензора с точностью 0,2° (по оценке дисперсии). Базовый алгоритм не справился с данной задачей, тогда как разработанный алгоритм обеспечил желаемую точность.

Оптимизировав параметры робастной статистики для нахождения теллурического тензора и включив в обработку магнитные каналы, удалось повысить качество получаемых частотных зависимостей.

Поправки за смещение оценок. В случае обработки без удаленной базы в описываемом алгоритме используется процедура, основанная на идеях, предложенных в работе А. Мюллера [Ми11ег, 2000]. При вводе поправок делается ряд допущений. Предполагается, что дополнительные компоненты тензора малы по сравнению с главными и существует линейная связь между величиной смещения оценки и уровнем шума в магнитных каналах:

-'xyi

-7" а

^ху ^ху^ху

\-СоЩ{Ну,Щ\ 1 — CoAQ (Нх , Ну )г

(4)

где 2^ — смещенная оценка модуля импеданса для временного интервала I, 2^, — истинное значение модуля импеданса, аху — линейный коэффициент, СоЩ(Ну,Н!у\ — когерентность между компонентой Ну и ее предсказанным значением (восстановленным по электрическому полю), СоИ$(Нх,Ну)( — когерентность между компонентами Нх и Н... х у

В этом уравнении два неизвестных — 2]^ и аху. Для их нахождения для разных временных интервалов выполняется оценивание с использованием робастных процедур, описанных выше, и рассчитываются когерентности для полученного решения. Затем составляется переопределенная система уравнений, из которой определяются величины смещений. Как показало тестирование, эту процедуру можно успешно использовать для

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

Использование всех каналов удаленной базы. Как известно, оценки рассматриваемых передаточных функций являются смещенными из-за неизбежного присутствия шума в наблюденных полях [Goubau, 1978]. Довольно простое решение этой проблемы было предложено Т.Д. Гэмблом (1978) в США и И.А. Безруком и В.О. Лахти-новым в СССР (1979) [Жданов, 1986]. Влияние шумов можно уменьшить, вводя в обработку поле, измеренное на некотором расстоянии от точки наблюдения (на удаленной базе). Такой способ борьбы с помехами хорошо себя зарекомендовал [Jones, 1989]. Удаленную базу размещают так, чтобы обеспечить достаточную линейную связь между полезным сигналом на основной точке измерения и на удаленной базовой точке и одновременно избежать коррелированности шумов между этими точками.

Решение для тензора импеданса с использованием удаленной базы имеет следующий вид [Gamble, 1979]:

yZyx

■"уу У

ExRx EyRx

ExRy EyR*yJ

HxRx HyRx

HxRy

HyRyj

(5)

где Rx и Ry — горизонтальные компоненты магнитного поля на удаленной базовой точке.

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

В рамках метода удаленной базы часто одновременно с магнитным полем измеряются и компоненты электрического поля. Однако в боль-

шинстве случаев оказывается, что электрические каналы зашумлены сильнее магнитных, и их использование в рамках общепринятых схем решения с удаленной базой, аналогичных (5), приводит к еще большему смещению оценок.

Предложенный автором статьи алгоритм, названный Remote E+H, позволяет использовать удаленные электрические каналы и не приводит при этом к смещению оценок.

Пусть есть M спектральных оценок каждой компоненты поля; Exi, Eyi, Hx,, Hyi — магнитные и электрические компоненты в точке наблюдения;

Rxv E

базы 77 базы

E° — магнитные и электрические

компоненты в базовой точке и .

Сначала выполняется классическая обработка с удаленной базой:

гЛ-1

7 7

"XX ^ху

\Zyx

"УУ )

ExRx ExRy

F R F R

HXRX HxRy

yHyRx HyRy j

где — компоненты тензора импеданса, а величины вида АВ* — авто- и кросс-спектры, т.е. попарные произведения спектральных оценок, осредненные с использованием весов, полученных с применением робастных процедур:

_ м

. (6)

1=1

1000-

о

о К X

U

&

©

-120

1000 10 000

Рис. 4. Сравнение результатов робастной обработки, полученных с использованием программы 88МТ2000 (а)

и рассматриваемого алгоритма (б)

На следующем этапе используются электрические каналы удаленной точки для подавления шума в электрических каналах локальной точки. Для этого сначала находится устойчивая робастная оценка теллурического тензора, определяющего связь между электрическим полем в локальной и базовой точках:

ßJlOK _ ^базы

^лок _ ^базЫ'

rp I г базы у

XX "г" "у *ху

гт< I р базы у ' *ух ' ¿-'у *уу

(7)

Тензор T вычисляется аналогично Z, как

X

(8)

влению путем задания малых значений весов wf, отвечающих зашумленным реализациям записи.

Тестирование алгоритма. Проверка эффективности описываемого алгоритма выполнялась на наборе МТ-данных разного качества. Сравнение проводилось в основном с результатами, получаемыми с помощью программного обеспечения «Phoenix Geophysics» SSMT 2000.

В случае, когда первичные записи МТ-поля характеризуются достаточно высоким соотношением сигнал/шум, результаты, получаемые в режимах простой одноточечной обработки и обработки с использованием магнитных каналов удаленной базы с помощью алгоритма SSMT 2000 и рассматриваемого алгоритма, оказываются очень близки. В то же время в условиях низкого соотношения сигнал/шум предложенный алгоритм часто оказывается эффективнее, чем алгоритм SSMT 2000.

Затем для каждой г'-й оценки спектров рас считываются невязки

рк, Ом ■ м 1000

_ рлок _ рбазыт" _ рбазы у rxi 1 хх f-'yi 1 xy 5

_ рлок _ рбазыт^ _ трбазыгр ryi f-'yi f-'xi 1 ух f-'yi 1 yy •

(9)

Затем происходит расчет веса каждой оценки по схеме, предложенной П. Хьюбером [Huber, 1981]. Веса задаются обратно пропорциональными невязкам:

1 И < г0

и

(10)

где |r | — модуль невязки, r0 — робастная оценка масштаба.

Наконец, авто- и кросс-спектры в уравнении (5) пересчитываются с учетом новых весов _ м

A£*=Zw?vt4B*, (11) ¡=i

и по ним находятся новые оценки импеданса. Эта процедура повторяется несколько раз до сходимости, при этом с каждой итерацией происходит уточнение весов ytf и vif.

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

1000 10 000

рк, Ом ■ м 1000

mr

1000

10 000

Рис. 5. Сравнение результатов, полученных с использованием трех различных режимов обработки: 1 — стандартная робастная обработка (Single); 2 — обработка с использованием магнитных каналов удаленной базы (Remote H); 3 — обработка с использованием 4-х каналов удаленной базы (Remote H+E)

В качестве примера сопоставлены результаты обработки МТ-данных с низким соотношением сигнал/шум в диапазоне периодов выше 1 с. На рис. 4 видно, что результаты, полученные с использованием описываемой системы обработки, характеризуются значительно меньшим значением дисперсии реализаций и меньшим уровнем их смещения.

Эффективность алгоритма Remote E+H показана на рис. 5, где приведены результаты разных типов обработки для записи с высоким уровнем шума в электрических каналах. Обработка Remote E+H в данном случае позволяет повысить точность получаемых кривых магнитотеллурического зондирования, что выражается в уменьшении доверительных интервалов и более точном выполнении дисперсионных соотношений.

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

Жданов М.С. Электроразведка М.: Недра, 1986. 316 с.

Семенов В.Ю. Обработка данных магнитотеллури-ческого зондирования. М.: Недра, 1985.

Chave A.D., Thomson D.J., Ander M.E. On the robust estimation of power spectra, coherences, and transfer functions // J. Geophys. Res. 1987. N 92. P. 633-648.

Egbert G.D., Booker J.R. Robust estimation of geomagnetic transfer functions // Geophys. J. R. Astr. Soc. 1986. N 87. P. 173-194.

Egbert G.D., Livelybrooks D. Single station magneto-telluric impedance estimation: coherency weighting and the regression M-estimate // Geophisics. 1996. N 61. P. 964-970.

Gamble T.D., Goubau W.M., Clarke J. Magnetotellu-rics with a remote magnetic reference // Geophysics. 1979. N 44. P. 53-68.

Goubau W.M., Gamble T.D., Clarke J. Magnetotelluric data analysis: removal of bias // Geophisics. 1978. N 43. P. 1157-1166.

Заключение. Разработан алгоритм (система) обработки магнитотеллурических данных, обладающий высокой устойчивостью по отношению к интенсивным электромагнитным шумам, присутствующим на записях МТ-поля. Особенность алгоритма — применение специальных схем мультипликативного взвешивания данных при вычислении разных передаточных функций, а также возможность одновременного использования как электрических, так и магнитных каналов удаленной базовой точки. В алгоритме используются различные методы борьбы со смещением спектральных оценок, включая робастную М-оценку, метод «складного ножа», использование удаленной базы, ввод поправок за смещение спектральных оценок. Предложенный алгоритм демонстрирует высокую эффективность при обработке МТ-данных с низким соотношением сигнал/шум.

Jones A.G., Chave A.D., Egbert G et al. A comparison of techniques for magnetotelluric response function estimation // J. Geophys. Res. 1989. N 94. P. 14 201-14 213.

Huber P.J. Robust Statistics. New York: John Wiley & Sons, 1981.

Larsen J.C. Transfer functions: Smooth robust estimates by least squares and remote reference methods // Geophys. J. 1989. N99. P. 655-663

Muller A. A new method to compensate for bias in magnetotellurics // Geophys. J. 2000. N 142. P. 257-269.

Sims, W.E., Bostick F.X., Smith H.W. The estimation of magnetotelluric impedance tensor elements from measured data // Geophysics. 1971. N 36. P. 938-942.

Zorin N, Epishkin D. Yakovlev A. A telluric method for natural field induced polarization studies // J. Appl. Geophis. 2015. N 23.

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