CC BY
58
УДК 520.628, 517.9
О. М. Лебедь, А. В. Ларченко, С. В. Пильгаев, Ю. В. Федоренко
РЕШЕНИЕ УРАВНЕНИЯ ТИПА СВЕРТКИ С ПОМОЩЬЮ ИНВЕРСНОГО ФИЛЬТРА
Аннотация
Приведен метод решения уравнения типа свертки, позволяющий с помощью цифрового инверсного фильтра преобразовать отсчеты АЦП выходных напряжений антенных усилителей электрической и магнитных компонент поля в их значения в А/м и В/м в ограниченной полосе частот.
Ключевые слова:
компоненты электромагнитного поля, метод, инверсный фильтр.
O. M. Lebed, A. V. Larchenko, S. V. Pilgaev, Yu. V. Fedorenko
DECONVOLUTION BY AN INVERSE FILTER
Abstract
The deconvolution method is discussed. It uses an inverse filter to transform ADC samples of bandlimited output voltage of sensors to electric and magnetic field values in A/m and V/m units.
Key words:
components of the electromagnetic field, method, inverse filter.
Введение
При обработке данных цифровой регистрации компонент электромагнитного поля, например при анализе электромагнитных возмущений, создаваемых удаленными молниевыми разрядами, сталкиваются с тем, что отсчеты аналого-цифрового преобразователя (АЦП) на выходе регистрирующей системы являются сверткой напряженности электромагнитного поля и импульсной функции регистратора. При анализе компонент поля возникает обратная задача восстановления их значений по отсчетам данных, записанных регистратором.
Операция, обратная свертке, не является устойчивой, поскольку и импульсная функция регистратора, и экспериментальные отчеты данных в принципе неточны. При решении любой обратной задачи, в том числе и интегрального уравнения типа свертки, возникают три основных вопроса: 1) существует ли решение; 2) если решение существует, то является ли оно единственным; 3) устойчиво ли решение, т. е. приводят ли малые изменения исходных данных к малым изменениям решения. Если решение существует и оно единственное и устойчивое, то задача называется корректно поставленной. В противном случае, задача называется некорректно поставленной или некорректной. Задачи восстановления геофизических сигналов, как правило, являются некорректно поставленными, так как не выполняется ни одно из вышеперечисленных условий.
Для решения некорректных задач А. Н. Тихонов предложил универсальный метод, получивший в отечественной и зарубежной литературе название «метод регуляризации Тихонова» [1]. Регуляризация решения состоит в построении семейства обратных операторов, зависящих от некоторого
101
числового параметра а (параметра регуляризации). Этот метод, безусловно, заслуживает внимания, поскольку он не использует априорной информации и требует лишь выбора параметра регуляризации. Однако при обработке длинной, в пределе бесконечной, последовательности данных, а также при работе в реальном времени, что является спецификой геофизических измерений, данный метод становится слишком сложно реализуемым и требующим большого количества времени для поиска решения.
Следует отметить, что задача обработки данных компонент поля обладает рядом особенностей. Во-первых, функции передачи аналоговой части регистраторов электрической и магнитных компонент являются дробнорациональными функциями. Это очень важно, так как для такой функции передачи можно построить фильтр в z-области, позволяющий вести обработку потоков данных. Во-вторых, анализ временных зависимостей компонент поля ведется в ограниченной полосе частот, ширина которой меньше полосы частот регистратора.
В данной работе приведен предложенный авторами метод, позволяющий с помощью цифрового инверсного фильтра преобразовать отсчеты АЦП выходных напряжений антенных усилителей электрической и магнитных компонент поля в их значения в А/м и В/м в ограниченной полосе частот. Проведен анализ ошибок преобразования, причем в качестве эталонного сигнала использован импульс электромагнитного возмущения, создаваемого молниевым разрядом вблизи экватора.
Инверсный фильтр для дробно-рациональной функции передачи
Известно, что функция передачи регистраторов электрической и магнитных компонент по определению является дробно-рациональной функцией вида P(s)IQ(s), где P и Q - полиномы, а s - комплексная переменная. Для дробно-рациональной функции свойственно сохранение каузальности при переходе из s- в z-область. Еще одним ее свойством является то, что основные ошибки в определении реальной передаточной функции сосредоточены на ее краях. В этом случае, если взять сигнал в полосе, края которой будут достаточно далеко отстоять от краев передаточной характеристики регистрирующей системы, влияние ошибок, возникающих за счет неточности определения передаточной характеристики на краях, можно будет свести к минимуму даже без применения метода регуляризации.
Функциональная схема предлагаемого инверсного фильтра приведена на рис.1. Здесь Ej - измеряемая компонента поля; i=x ,y, z; Pj(s)IQj(s) - функции передачи каналов регистратора в виде дробно-рациональной функции аргумента; s=jm; U(s)IV(s) - унифицированная функция передачи; Wjnv(s) -инверсный фильтр для i-й компоненты.
Регистратор
Ei(s) ЫЗ) г
Компоненты " Qi(s) Компоненты
поля в А/м и В/м поля в ед. АЦП ^
Инверсный фильтр
U(S)p(n\
v(s)b*>
Компоненты поля в А/м и В/м
Рис.1. Схема метода решения обратной задачи восстановления компонент поля
102
Суть метода заключается в том, что сначала выбирается одинаковая для всех компонент поля, регистрируемых на всех станциях, передаточная функция H(s) вида H(s)=U(s)/V(s), s=lKjf=ja, где f- частота; U(s) и V(s) - рациональные функции s. Число полюсов у нее должно быть как минимум на два больше, чем у передаточной функции аналоговой части системы сбора. Два добавочных полюса нужны для ограничения полосы частот выходных сигналов. Выбором этой передаточной функции определяется полоса частот, в которой будет вестись последующий анализ измерений компонент поля. По определенным ранее передаточным функциям измерительных каналов компонент поля Fx(s)=Px(s)/Qx(s), Fy(s)=Py(s)/Qy(s) и Fz(s)=Pz(s)/Qz(s) далее рассчитываются передаточные функции инверсных фильтров каждого канала:
W inv
r x, у, z
(s)
H(s) _ U(s) Qx,y, z (s)
Fx, у z (s) V(s)' Px, y z (s)
Функция передачи инверсного фильтра Wx y, zmv(s) - это дробнорациональная функция переменной s, число нулей которой как минимум на два меньше, чем число полюсов. Как известно [2], полюса и нули такой функции можно преобразовать в коэффициенты БИХ-фильтра Wx y, zinv(z) при помощи, например, билинейного или любого другого известного преобразования из s-области в z-область. Полученный БИХ-фильтр удобен для обработки длинных последовательностей данных. В принципе, такой фильтр даже может быть включен на входе записывающего устройства для преобразования отсчетов АЦП в отсчеты компонент поля в реальном масштабе времени. После обработки сигналов компонент поля инверсным фильтром результирующие амплитудночастотные и фазо-частотные характеристики H(s) компонент одинаковы у всех измерительных каналов на всех станциях, а значения отсчетов сигналов компонент поля представляют измеряемые величины Hx, Hy и Ez в А/м и в В/м.
Выбор метода синтеза и апробация цифрового фильтра
Преобразование инверсного фильтра из s-области в z-область представляет большой интерес для ускорения обработки сигналов и особенно для обработки потоков данных в системах реального времени. Существует три наиболее часто применяемых метода синтеза цифровых фильтров (преобразования из s- в z-область), каждый из которых обладает своим рядом достоинств и недостатков [2]. В их число входят так называемый метод выборки-хранения (zero order hold - ZOH), билинейное преобразование и метод согласованного z-преобразования (matched pole/zero method -MPZ). Поскольку все эти методы могут быть применены для синтеза цифровых фильтров из передаточных характеристик измерительных каналов, перед авторами стояла задача выбора оптимального метода из перечисленных выше.
На рисунке 2a-c приведены АЧХ, ФЧХ и время групповой задержки, вносимой стандартной передаточной характеристикой в s-области и эти же характеристики, полученные тремя методами в z-области. Из рисунка следует, что ни один из методов преобразования инверсного фильтра из s- в z-область не показал хороших результатов. Особенно большие погрешности возникают вблизи частоты Найквиста (здесь она равна 256 Гц). Можно утверждать, что наблюдаемые отклонения этих характеристик от характеристик аналогового фильтра настолько велики, что не позволят с достаточной степенью точности восстановить значения компонент поля на входе измерительной системы.
103
Рис.2. Сравнение методов синтеза цифрового фильтра в частотной и временной области:
a - АЧХ; b - ФЧХ; c - время групповой задержки; d - модельный импульс. Вертикальные линии - выделенная инверсным фильтром полоса частот
Поскольку данная работа основана на исследовании импульсных сигналов электромагнитных возмущений во временной области, сравним модельный импульс, полученный при прохождении сигнала атмосферика через фильтр со стандартной характеристикой H(s), с сигналом, пропущенным сначала через аналоговую часть регистратора c функцией передачи F(s), а затем через инверсный БИХ-фильтр, полученный из Wnv(s) с помощью перечисленных выше трех методов перехода из s- в z-область. Результат сравнения представлен на рис.2d. Здесь наглядно показано, какие искажения будут присутствовать в восстановленном сигнале, если воспользоваться любым из обсуждаемых методов перехода из s- в z-область.
Известно, что цифровые фильтры вносят наибольшие отклонения от аналогового фильтра-прототипа в полосе частот, верхняя граница которой расположена вблизи частоты Найквиста Fq=Fs/2 (см. рис.2). Очевидно, что если бы частота дискретизации была настолько высока, что частота Найквиста далеко отстояла от верхней частоты спектра данных, цифровой инверсный фильтр был бы близок по свойствам к аналоговому и результат его работы был бы практически неотличим от результата работы аналогового фильтра. Отсюда следует, что эти отклонения могут быть устранены повышением частоты
104
дискретизации перед применением инверсного фильтра, преобразованного в БИХ-фильтр.
Для того чтобы выбрать, во сколько раз надо повысить частоту дискретизации, оценим ошибку, вносимую каждым из трех методов преобразования сигнала из s- в z-область на примере восстановления модельного импульса, последовательно увеличивая частоту дискретизации. Мы повысили частоту дискретизации сигнала в 2, 4, 8 и 10 раз и для каждого из методов рассчитали относительную ошибку восстановления импульсного сигнала. Результаты расчетов сведены в таблицу.
Таблица
Относительная ошибка восстановления импульсного сигнала при преобразовании из s- в z-область для различных частот дискретизации, %
Метод преобразования/коэффициент увеличения частоты дискретизации 1 2 4 8 10
ZOH 90 27 9.5 5.9 4.9
Билинейное преобразование 53 11 2.5 0.53 0.4
MPZ 20 21 15 9.9 8.8
Из общих соображений ясно, что чем меньше частота дискретизации, тем более «плотную» запись можно создать, тем больше данных можно записать на один диск и тем менее мощный компьютер потребуется для обработки такой записи. Кроме того, чем ниже частота дискретизации, тем больше искажений возникнет при обработке в силу свойств цифровых фильтров. Из таблицы видно, что при увеличении частоты дискретизации в 2 раза ни одним из методов не достигается ошибки даже ниже 10%. Увеличив частоту дискретизации в 4 раза, мы наблюдаем хорошее совпадение модельного сигнала и сигнала, восстановленного с помощью билинейного преобразования. Здесь ошибка составляет всего 2.5%. Дальнейшее увеличение частоты дискретизации приводит к избыточно низкой ошибке для билинейного преобразования и к допустимой ошибке для двух оставшихся методов. Таким образом, можно сделать вывод, что оптимальным решением задачи преобразования инверсного фильтра Wm’(s) из s- в z-область для импульсных сигналов является использование билинейного преобразования и повышение частоты дискретизации в 4 раза.
В заключение, на рис.3 (а-с) показаны АЧХ, ФЧХ и время групповой задержки инверсного фильтра и БИХ-фильтра, полученного преобразованием инверсного фильтра из s- в z-область с помощью билинейного преобразования, соответственно. На рис.3 (d) показан модельный сигнал и сигнал, восстановленный с помощью указанного БИХ-фильтра. Здесь частота дискретизации увеличена в 4 раза. Операция повышения частоты дискретизации выполнялась стандартным способом, описанным, например, в [2].
Из рисунка видно, что выбранный метод преобразования инверсного фильтра из s- в z-область и выбранный коэффициент повышения частоты дискретизации позволили получить хорошее совпадение результатов обработки сигнала атмосферика цифровым фильтром с результатами обработки того же сигнала аналоговым прототипом цифрового фильтра.
105
Рис.3. Сравнение инверсного фильтра и БИХ-фильтра, полученного с помощью билинейного преобразования: a - АЧХ; b - ФЧХ; c - время групповой задержки; d - модельный импульс. Вертикальные линии - выделенная инверсным фильтром полоса частот
Выводы
Авторами предложен метод решения обратной задачи восстановления физической величины из отсчетов аналого-цифрового преобразователя выходного напряжения геофизического датчика, функция передачи которого является дробно-рациональной функцией. Описан способ расчета цифрового инверсного фильтра, который позволяет вести обработку как непрерывного потока данных в режиме реального времени, так и выделенных сегментов данных. Показано, что оптимальным преобразованием из s- в z-область является билинейное преобразование. Произведены оценки ошибок перехода из s-в z-область. Предложенный инверсный фильтр применен для анализа электромагнитных возмущений, создаваемых удаленными молниевыми разрядами.
Литература
1. Тихонов А. Н. Методы решения некорректных задач. М.: Наука, 1986. 288 с.
2. Сергиенко А. Б. Цифровая обработка сигналов. СПб.: Питер, 2003. 604 с.
Сведения об авторах Лебедь Ольга Михайловна,
младший научный сотрудник, Полярный геофизический институт, г. Апатиты, olgamihsh@yandex.ru
Ларченко Алексей Викторович,
младший научный сотрудник, Полярный геофизический институт, г. Апатиты, alexey.larchenko@gmail.com
Пильгаев Сергей Васильевич,
младший научный сотрудник, Полярный геофизический институт, г. Апатиты. pilgaev@pgia.ru
Федоренко Юрий Валентинович,
к.физ.-мат.н., доцент, заведующий сектором, Полярный геофизический институт, г. Апатиты, yury.fedorenko@gmail.com
106
