Эффективные алгоритмы цифровой обработки сигналов в задаче магнитотеллурического зондирования

Помехоустойчивость инвариантной системы передачи информации при наличии слабых корреляционных связей и собственных шумов генераторного оборудования. — Омский научный вестник, выпуск. — 2008. — № 3 (70). — С. 122-126.

6. Левин, Б.Р. Теоретические основы статистической радиотехники. — 3-є изд. - М. : Радио и связь, 1989. — 656 с.

7. Теплов, Н.Л. Помехоустойчивость систем передачи дискретной информации. — М. : Связь, 1964. — 359 с.

АЛГАЗИН Евгений Игоревич, кандидат технических наук, старший преподаватель кафедры общей электротехники Новосибирского государственного технического университета.

КОВАЛЕВСКИЙ Артем Павлович, кандидат физико-

математических наук, доцент кафедры высшей математики Новосибирского государственного технического университета.

МАЛИНКИН Виталий Борисович, доктор технических наук, профессор кафедры многоканальной электросвязи и оптических систем Сибирского государственного университета телекоммуникаций и информатики. КАСАТКИНА Елена Геннадьевна, кандидат технических наук, доцент кафедры общей электротехники Новосибирского государственного технического университета.

630092, г. Новосибирск, ул. К. Маркса 20, к. 536

Дата поступления статьи в редакцию: 19.02.2009 г.

© Алгазин Е.И., Ковалевский А.П., Малинкин В.Б., Касаткина Е.Г.

удк519 254 А. А. ЛАВРУХИН

К. В. ЛОБОВ

Омский государственный университет путей сообщения Омский государственный технический университет

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

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

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

Метод магнитотеллурического зондирования (МТЗ) земной коры заключается в регистрации, обработке и интерпретации естественных электромагнитных полей на поверхности Земли с целью получения информации о свойствах и характере разрезов [1]. Получение полной картины в геологоразведке требует проведения измерений на местности в узлах сетки, расположенных на расстоянии 200-500 метров. Регистрируемые поля имеют частоты от 10"4 Гц, поэтому каждое такое измерение длится от нескольких часов до суток и временные затраты на измерения очень велики. Другой нерешенной задачей является практическое отсутствие возможности в полевых условиях оценивать качество получаемых данных, чтобы определять их необходимый объем, достаточный для достоверной интерпретации. Обычно вся обработка данных ведется в лабораторных условиях после проведения полевых работ.

Таким образом, стоит задача сокращения расхо-

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

В обратной задаче МТЗ [ 1 ] по зарегистрированным ортогональным сигналам напряженностей электрических Ех, Еу и магнитных Нх, Ну компонент поля требуется получить параметры геоэлектрической модели. Задача состоит из двух частей. В первой по

ОМСКИЙ НАУЧНЫЙ ВЕСТНИК Ж 2 (80). 200» РАДИОТЕХНИКА И СВЯЗЬ

РАДИОТЕХНИКА И СВЯЗЬ ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 2 (80), 2009

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

В соответствии с моделью Тихонова-Каньяра [2] и классической постановкой магнитотеллурической задачи, компоненты искомого тензора импеданса 1ХХ, 2ху, 2ух, связаны с мгновенными комплексными значениями компонент поля Ех, Еу, Нх, Ну линейной системой уравнений. Поскольку необходимо получить частотные зависимости компонент импеданса, система уравнений записывается для значений фиксированных частотных составляющих сигналов.

Определение тензора импеданса состоит из нескольких этапов, первые два из которых рассматриваются в статье:

— предварительной фильтрации зарегистрированных сигналов;

— разбиения сигналов на частотные составляющие;

— отбраковки слабокоррелирующих и несинфазных временных участков;

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

Качество оценок импедансных кривых можно повысить, если на первом этапе обработки по возможности исключить известные в сигналах помехи. Искажение сигнала происходит в результате наложения на него импульсных помех техногенного характера типа «выброс» и «ступенька». Для первоготипапомехи один или несколько дискретных отсчетов смещаются относительно тренда сигнала, а для второго - резко изменяется постоянная составляющая сигнала.

Для устранения помех типа «выброс» используется метод, основанный на предсказании по авторегрессионной модели. Он был разработан по аналогии с методом, предложенным в статье [3], но отличается большей простотой и вычислительной эффективностью. Дискретный сигнал х(л) разбивается на последовательно расположенные сегменты равной длины по Л/9дискретных отсчетов. Каждому сегменту с номером т ставится в соответствие среднее значение сигнала хш и среднеквадратичное отклонение от, которые определяются по трем соседним сегментам с номерами (т - 1), т и (т + 1). Если для значения сигнала х(л) из т-го сегмента выполняется условие

|х(п)-

где а — положительная константа, то считается, что п-й отсчет сигнала искажается помехой, и его следует заменить другим значением, предсказанным по предыдущим отсчетам.

Сигнал предлагается описывать авторегрессионной моделью, поэтому дискретное значение сигнала х(л) определяется выражением

х

~1

а„(*)-х(л-А) + и(л),

по предыдущим р отсчетам, а также невязка

е

между действительным и предсказанным значениями.

Подстановка выражений (2) и (3) в формулу (4) приводит к равенству е(л) = и (л). Поскольку сигнал и (л) в выражении (2) является белым шумом, то и ошибка предсказания е(л) — белый шум. Таким образом, если определить коэффициенты авторегрессии ар(к) по предыдущим N значениям сигнала с использованием одного из известных алгоритмов, то для устранения помехи отсчет х(л) следует заменить на х(л).

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

Экспериментальные исследования с использованием полевых данных показали, что оптимальные значения константы а находятся в пределах от 2,5 до 3, а наилучшим способом определения коэффициентов авторегрессии является метод Берга, дающий срав-нимые с ковариационным и автоковариационным методами результаты, но наиболее эффективный по скорости вычисления. Порядок р рекомендуется брать из диапазона от 40 до 100, а массив данных для определения коэффициентов авторегрессии должен иметь длину N от 100 до 200.

Для идентификации и удаления помех типа «ступенька» предлагается использовать следующий метод Сигнал х(0 пропускается через фильтр скользящего среднего (путем вычисления среднего значения внутри скользящего окна шириной 2ГСС), в результате чего получается низкочастотная составляющая у((). Особенность такой фильтрации заключается в том, что в окрестности помехи типа «ступенька» полученная функция времени у(() имеет измененный наклон, который начинается в момент времени <п — Гсс и заканчивается во время + Гсс, где (и — момент действия помехи в сигнале х(0- На отрезке времени [/п — Гсс; £„ + Гсс] наклон у(0 изменяется незначительно. Производная по времени у(() в моменты <п — Гсс и + Тсс скачкообразно меняется. Вторая производная у(1) в обозначенные моменты времени имеет резкие выбросы с противоположными знаками. Поскольку при определении производных помехи в сигнале усиливаются, предлагается использовать сглаживание по 1 отсчетам в каждый момент численного дифференцирования. Тогда для л-го дискретного момента времени первая и вторая производные будут определяться соответственно по формулам

где ар(к) — коэффициенты авторегрессии; и (л) -случайная составляющая модели, вносящая нестаци-онарность в наблюдаемый процесс.

Уравнению (2) соответствует формула для предсказания значения сигнала

х

У

ЧА-1 к-1 )/

У

Ч*-1 *-1 //

где Г0 — величина шага дискретизации.

Зависимости от времени скользящего среднего у(0 полевого сигнала х(0 и ее производных у(<) и у(<) показаны на рис. 2. Исправленный сигнал с исключенной шумовой составляющей обозначен х(л).

Положение помехи типа «ступенька» можно определить по описанной выше методике, используемой для нахождения помех типа «выброс», так же пугем разбиения сигнала у(/) на сегменты, определения для каждого из них среднего значения и средне-

квадратичного отклонения. Если найден локальный экстремум функции у(<), превосходящий по модулю некоторое пороговое значение, и на расстоянии 2Та от него существует второй локальный экстремум, сравнимый по модулю, но противоположного знака, тогда можно сделать вывод о том, что в этом месте присутствует помеха.

Чтобы определить величину сдвига, можно интерполировать значения функции у(<) слева и справа от помехи по ее производным у(£) за пределами отрезка К. - Тсс'1„ + гсс]. Тогда сдвиг Ду вследствие действия помехи будет определяться по формуле

Лу = У (<2) - У (* 1)'- У (*2) • Гсс - У(<|)  Гсс,

где 1Ш тсс Т0, <в + Тсс + Т0.

Устранение помехи производится сдвигом следующего, начиная с момента времени 7^, сигналах(/) на величину Ду.

Исследования показали, что хорошие результаты идентификации помехи получаются при значении параметра Гсс скользящего окна из диапазона от 4 до 100 дискретных отсчетов и параметре сглаживания 1, равном от 2 до 8.

На следующем этапе обработки данных МТЗ необходимо представить сигналы в частотной области, либо выделить их частотные составляющие как функции времени. Это может быть сделано на основе преобразования Фурье периодограммным или кор-релограммным способом [4], вейвлет-анализа [5], или узкополосной фильтрации с последующим определением амплитудных и фазовых характеристик на дискретном наборе частот [6].

Исходя из специфики задачи МТЗ, частотные характеристики (ЧХ) зарегистрированных сигналов нужно получать для широкого диапазона частот, обычно от 0,001 до 10 Гц, а иногда и шире, с логарифмическим распределением по частотам. Рассматриваемые сигналы далеки от стационарных, и их параметры (амплитуда, фаза) на фиксированных частотах могут значительно изменяться в течение нескольких периодов. Методы, основанные на Фурье-анализе, не подходят для спектральной обработки таких сигналов, поскольку получение широкой ЧХ требует преобразования Фурье для большого массива данных, что предполагает некоторую стационарность характеристик рассматриваемого процесса.

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

арифметических операций над комплексными числами, где I — число масштабов преобразования (соответствует количеству фиксированных частот); N — число дискретных отсчетов сигнала.

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

Рис. 4. Алгоритм первичной обработки данных

,« ——— Нпашумлсннмй сигнал

- — - Зашумленный сигнал

-----Сигнал с исключенной

помехой

20 21 22 23 24 25 26

Рис. 1. Исключение помехи типа «выброс»

I

X

У

У

У

х

I

Рис. 2. Исключение помехи типа «ступенька»

4000 4200 4400 ^^ 4600

Рис. 3. Получение мгновенных амплитуды и фазы сигнала в узкой полосе частот

РАДИОТЕХНИКА И СВЯЗЬ ОМСКИЙ НАУЧНЫЙ ВЕСТНИК № 2 (80). 2009

к

имеет преимущество по скорости вычисления перед предыдущим подходом. В его основу положены два следующих алгоритма.

Первый заключается в преобразовании сигнала путем сдвига его частотной характеристики, как это делается при амплитудно-фазовой модуляции в радиосвязи [7]. Допустим, в дискретном сигнале х(л) требуется выделить полосу частот шириной 2Д/т в окрестности частоты /т. Рассмотрим вспомогательный сигнал у(л), получаемый путем умножения всех отсчетов исходного сигнала на комплексную экспоненту:

v(n) = х(п)ехр(—j-2nnfm/f0),

где л — номер дискретного отсчета:/0 - частота дискретизации; / — мнимая единица.

Сигнал у(п) является комплексным, и поэтому имеет несимметричную спектральную характеристику. Выделение узкополосной составляющей хпф(л) с полосой [/т — Д/т; !т + Д/т] из сигнала х(л) будет эквивалентно выделению низкочастотной составляющей у||Ф(л) с частотами [ — Д/я1; Д/т]. Поскольку при этом симметричный диапазон частот [ — fm — Д/т; — (т + Д/т] исходного сигнала будет сдвинут в окрестность частоты - 2Д/т, то пропускание у(л) через низкочастотный фильтр оставит только основной лепесток спектра, и поэтому уменьшит амплитуду в два раза. Отфильтрованный сигнал 2-у||Ф(л) можно рассматривать, как модулирующий для единичной синусоиды частотой /т, дающий в результате такой модуляции отфильтрованный хПФ(л). Мгновенная амплитуда узкополосного хпф(л) определяется по моду-

лю сигнала унф в л-й момент времени:

а фаза - по формуле

ф1хпф(п)1 = argv1№(n) + nf„,/f0.

Второй алгоритм, который многократно используется в процессе получения частотных характеристик, осуществляет понижение частоты дискретизации (прореживание) в два раза. Для этого проводится фильтрация сигнала фильтром нижних частот и собственно прореживание. Используемый фильтр является нерекурсивным и не искажает фазовую характеристику сигнала. Частота пропускания фильтра равна/0/16, а частота задерживания — /„/^Подавление в полосе задерживания составляет 80 дБ. Порядок <7 оптимального фильтра с равнозначными пульсациями в полосах пропускания и задерживания, имеющего такие характеристики, равен 17. Если максимальная интересующая частота в сигнале равна /пи>, то признаком окончания прореживания будет попадание этой частоты в диапазон Г||мхе [Г0/16; Ги/8], поскольку, принимая во внимание критерий Найквис-та, частоту дискретизации рекомендуется принимать на порядок выше максимальной частоты сигнала (7).

Выполняется два этапа понижения частоты дискретизации. На первом этапе прореживанию подвергается исходный, очищенный от помех, действительный сигнал х(л). Признаком окончания будет выполнение условия Гт + Д^е (Г0/16; Г0/8]. После этого выполняется операция (9) сдвига спектра и получение комплексного у(л). На втором этапе прореживается сигнал у(л), пока не выполнится условие ДГгое [(0/16; Г0/8].

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

вать многократно — для выделения различных частот / в сигнале.

/71

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

_ Nqf. Ill ^

R. = -2- ! + _ + _ + _ + «/Vqr

1 2 I, 2 4 8 ) Ч-

Умножение на комплексную экспоненту займет не более чем

Я =n[i + - + - + - + ...|»2N { 2 4 8 )

операций над комплексными числами.

Следующий этап прореживания осуществляется отдельно для каждой частоты /т. Если считать, что количество фиксированных частот равно L и во всем диапазоне они распределены с логарифмическим шагом, тогда этот диапазон будет разбит на log2L отрезков, каждый из которых соответствует одной полученной после первого этапа частоте дискретизации. На каждом отрезке содержится £/1од21.фик-сированных значений частот. Для отрезка на старших частотах частоту дискретизации будем считать равной /0 исходного сигнала (наихудший случай, когда прореживание до умножения на комплексную экспоненту не проводилось). Тогда число арифметических операций, соответствующих этому отрезку, будет равным

NqL (, Ill ^ NqL

IL. =--— 1+ - + - + - + »——

™ 21og,U 2 4 8 ) log2L'

Соседнему отрезку более низких частот соответствует число операций

_ 1 NqL (, ,1 1 ,1 \ I NqL

Иг, —--------l*t—і—ь —к., я*---------

221og,U 2 4 8 ) 21og,L

Продолжая далее, получаем, что второй этап прореживания требует

о о о і , 1 1 1 ^ ATqrZ. 2 NqL

R,, = Л., + R.n + ...« 1-і—і—і—f*... -------я*------

™ ^ 1 2 4 8 Jlog2L log2L

арифметических операций над комплексными числами.

Таким образом, складывая полученные выражения для К,, Яе и Я2, получаем, что вычислительные затраты метода определяются формулой

Яф » N(q + 2) +

2NqL log2 L

Сравнивая выражения (8) и (12) с учетом того, что значение N обычно лежит в пределах 104 — 107,1 — от 100 до 10000, а <7 = 17, можно получить, что выиг-

рыш по количеству операций в случае использования полосовой фильтрации составляет от 2 до 9 раз. Это означает, что предлагаемый способ обработки будет работать в несколько раз быстрее, чем основанный на вейвлет-анализе.

Пример определения временных зависимостей мгновенных амплитуды и фазы узкополосного сигнала показан на рис. 3.

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

Приведенные методы выделения и устранения помех, а также процедуры расчета импедансных кривых, основанные на полосовой фильтрации, были опробованы в вычислительных экспериментах с использованием синтетических данных СОМОАТ [3]. Полученные результаты подтверждают пригодность методов в задачах оперативной оценки качества сигналов измерений МТЗ. Разработанные алгоритмы реализованы в виде программного обеспечения, апробацию которого в составе программно-аппаратного комплекса планируется провести при поиске полезных ископаемых в северных районах Омской области.

Библиографический список

1. Бердичевский М. Н. Магнитотеллурическое зондирование горизонтально-однородных сред / М.Н.Бердичевский, В.И.Дмитриев. — М. : Недра, 1992. — 249 с.

2. Электроразведка / под ред. В. К. Хмелевского,

И. Н. Модина, А. Г. Яковлева. — М. : Эксперт-диагностика, 2005. — 311 с.

3. Smirnov М. Yu. Magnetotelluric data processing with a robust statistical procedure having a high breakdown point // Geophysical Journal International, 2003. Vol. 152. P. 1-7.

4. Марпл-мл. С. Л. Цифровой спектральный анализ и его приложения / С. Л. Марпл-мл. — М. : Мир, 1990. — 584 с.

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

6. Оппенгейм А. Цифровая обработка сигналов / А Оппенгейм, Р. Шафер. — М.: Техносфера, 2007. - 856 с.

7. Скляр Б. Цифровая связь. Теоретические основы и практическое применение / Б. Скляр. — М. : Вильямс, 2004. — 1104 с.

ЛАВРУХИН Андрей Александрович, кандидат технических наук, старший преподаватель кафедры «Радиотехнические и управляющие системы» Омского государственного университета путей сообщения.

E-mail: lavruhinaa@gmail.com ЛОБОВ Константин Владимирович, ассистент кафедры «Средства связи и информационная безопасность» Омского государственного технического университета.

E-mail: yodobrom@gmail.com

Дата поступления статьи в редакцию: 04.05.2009 г.

© Лаврухин А.А., Лобов К.В.

УМ*21396 В. А. АРЖАНОВ

Омский государственный технический университет

СРАВНИТЕЛЬНАЯ ОЦЕНКА ЭФФЕКТА БЛОКИРОВАНИЯ РЕЗОНАНСНЫХ УСИЛИТЕЛЕЙ ПРЕСЕЛЕКТОРА

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

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

Прохождение сигнала через линейный тракт радиоприемного устройства (РГГУ) сопровождается различными нелинейными явлениями: нелинейные иска-жения модулирующей функции: изменение коэффициента передачи поддействием мешающего сигнала; перекрестные и интермодуляционные искажения: переход амплитудной модуляции в фазовую (амплитудно-фазовая конверсия) и фазовой в амплитудную (фазо- амплитудная конверсия).

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

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

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

т

= 1л -и:х = а, + а, • 1/„ + а2 • Щ + а3• и I

к-0

ОМСКИЙ НАУЧНЫЙ КСТНИК К* 2 (SO). 2009 РАДИОТЕХНИКА И СМЗЬ